- Email: [email protected]
Erato stability code
Computer Physics Communications 21(1981) 323 —371 © North-Holland Publishing Company
ERATO STABILITY CODE R. GRUBER, F. TROYON, D. BERGER
“,
L.C. BERNARD
~,
S. ROUSSET
~,
R. SCHREIBER
EURA TOM/Switzerland Fusion Association, Centre de Recherches en Physique des Plasmas, Ecole Polytechnique Fédérale de Lausanne, CH-1007 Lausanne, Switzerland
W. KERNER, W. SCHNEIDER Max-Planck-Institut für Plasmaphysik, Garching, EURA TOM-A ssociation, Fed. Rep. Germany
and K.V. ROBERTS Euratom/UKAEA Fusion Association, Cuiham Laboratory, Abingdon, Oxfordshire 0114 3DB, UK Received in final form 11 September 1980
PROGRAM SUMMARY Titleof program: ERATO
Keywords: plasma physics, ideal MHD, variational principle, finite elements, spectrum, instabilities, eigenfunctions
Catalogue number: ABVS Nature of the physical problem This computer code treats the stability of a Tokamak-like plasma described by the ideal linearized magnetohydrodynamic (MHD) equations [1]. The plasma is considered to be in an equilibrium state. After perturbation of such an equilibrium, the evolution of the normal modes can be calculated by linearizing the ideal MHD equations. Any mode (unstable or stable) of the spectrum can be examined.
Program obtainable from: CPC Program Library, Queen’s University of Belfast, Northern Ireland (see application form in this issue) Computer: CDC Cyber 170-720; Installation: Ecole Polytechnique Fédérale, Lausanne Operating system: NOS BE
Method of solution The variational form [2] of the 2D ideal MHD equations is treated by a finite hybrid element approach [3], which proves to be well suited to describe the features of the problem2Bx with sufficient accuracy 141. The eigenvalue problem Ax = w is solved by VEKIT (an inverse vector iteration), a subprogram of the block matrix library HYMNIABLOCK [5]. The code consists of 5 main programs (ERATO ito ERATO 5) linked together by disk files.
Programming language used: FORTRAN High speed store required: 35 000 words Number of bits in a word: 60 No. of magnetic tapes required: none Other peripherals used: line printer, disk files
Typical running time The running time is proportional to the number ofpoloidal flux surfaces Ni,, and to the cube of the number ofangular mtervals N~when N~>20. A typical case with 24 X 24 intervals in ‘I’ and x takes iS mm per eigenvalue on a CDC 6500.
No. of cards in combined program and test deck: 7901 Card punching code: EBCDIC
Present address: * IBM, Lausanne, Switzerland; ** General Atomic, San Diego, USA; *** Oak Ridge National Laboratory, USA; **** Zschokke SA Geneva. Switzerland. 323
Unusual features of the program The code uses subroutines from two program libraries, the
324
R. Gruber et al.
/ ERA TO stability code
utility routines of the CDC OLYMPUS package [6] and the eigenvalue package HYMNIABLOCK [5]. However, all the required subroutines are included in the ERATO package itself. ERATO is written in STANDARD FORTRAN [7], except for the use of the input facility NAMELIST which is available on most computers, and the CDC statement ENCODE which is used in ERATO 5. The plotting routines PLOTS, PLOT and SYMBOL used in ERATO S must be written to suit the local computer system and details oftheir arguments are given. ERATO also uses the CDC function subprogram DATE which picks up the actual date for plotting, and LOCF which gives the base address of an array. The programs ERATO 1—5 axe arranged so that they can readily be combined into a single program if sufficient storage is available. References [1] B.B. Kadomtsev, Reviews of plasma physics, vol. 2 (Consultants Bureau, New York, i966) p. 155. [2] W.A. Newcomb, Ann. Phys. 10 (1960) 232. [3] R. Gruber, J. Comput. Phys. 26 (1978) 379. [4] ~R.Gruber and F. Troyon, 3ème Colloque Intern. Méthodes de calcul scientifique et technique (Springer Verlag, Berlin, 1979); Lecture Notes Phys. 91(1979) 288. [5] R. Gruber, Comput. Phys. Commun. 20 (1980) 421. [6] M.H. Hughes, K.V. Roberts and G.G. Lister, Comput. Phys. Commun. 10 (1975) 167. [7] Standard FORTRAN Programming Manual, Computer Standard Series (National Computing Centre Ltd., Manchester, UK, 1970).
CONTENTS 1. Introduction 2. Equilibrium 2.1. General equations 2.2. Solovev equilibrium 2.3. Numerical equilibria 3. Stability problem 3.1. General axisymmetric problem 3.2. Choice of independent and dependent variables 3.3. Potential and kinetic energy in axisymmetric geometry 3.4. Vacuum energy in axisymmetric geometry 3.5. Additional n = 0 terms
4. Numerical approach 4.1. The marginal point of the continuous spectrum 4.2. The finite hybrid elements 4.3. Choice of the mesh 4.4. Spectral shift and its elimination 5. Organization and normalization 6. ERATO 1 6.1. Initialization 6.2. Choice of the mesh 6.3. Analytic treatment of the axis 6.4. Mapping 6.5. The normalized equilibrium quantities 6.6. Scaling law 6.7. Mercier’s criterion 6.8. Ballooning criterion 6.9. Definition of the betas 6.10. Surface fitting 7. ERATO 2 7.1. Treatment of the singularities 7.2. Discretisation 7.3. Vacuum energy matrix 7.4. Then = 0 terms 7.5. Resolution of the integrals 7.6. Definition of the shell 8. ERATO 3 8.1. Organization 8.2. Potential and kinetic energy 8.3. Symmetry conditions 8.4. Regularity and boundary conditions 9. ERATO 4 9.1. Eigenvalue problem 9.2. How to call VEKIT 10. ERATO 5 10.1. Diagnostics 10.2. Plotting 11. Test cases 11.1. Preset values (internal mode) 11.2. Test run output (kink mode) 12. Instructions to the user 12.1. Convergence properties 12.2. How to switch from equilibrium to ERATO 12.3. Choice of NAMELIST parameters 12.4. Control stream and data set 12.5. Dimensions 12.6. Space and time requirements
LONG WRITE-UP 1. Introduction Thermonuclear fusion research is profoundly concerned with the problem of macroscopic ideal MHD instabilities. When such instabilities occur they are generally very violent and global and lead to destruction of the confinement. The calculation of these instabilities has therefore become an important subject, and as a consequence of the great interest in Tokamaks, numerical stability codes have been written for axisymmetric toroidal geometry (2D).
R. Gruber etal.
/ ERA TO stability code
325
The first codes to be developed [1,2] solve the Newton equations of motion by finite differences. These codes work well for high pressure plasma but for weakly-growing modes the computer time becomes prohibitive and furthermore they are only able to find the most unstable mode. In the meantime two spectral codes (PEST from Princeton and ERATO) for a general axisymmetric equilibrium have been written and compared [3]. These two codes start from a variational form [4]. The existence of a continuous spectrum extending to the marginal point is very awkward for the numerical treatment. At this point, the mode is divergence-free and the B V operator vanishes. Weakly unstable modes, which are close to the marginal situation, are nearly incompressible and B V is very small. To overcome these problems, both groups use different finite elements for different components [5] to satisfy the incompressibility condition, and so prevent spectral pollution. The other condition (B V = 0) can be fulfilled by a Fourier analysis as implemented in PEST [6] or by a “finite hybrid element” approach [7] as used in ERATO. These hybrid elements have numerical advantages. They are piecewise constant in each mesh cell. The equilibrium quantities can also be chosen to be constant in each mesh cell; they only have to be given in the centre. The integration of the matrix elements is then trivial. The number of unknowns is smaller than for a regular finite element approach [8]. The convergence is generally from below, opposite to all Galerkin procedures. This enables us to well separate the unstable modes from the stable continuous spectrum and we never miss any unstable displacement. ERATO consists of 5 main programs (ERATO 1 to ERATO 5) linked together by disk ifies. For numerical equilibria, a sixth program, the equilibrium code, has to precede ERATO 1. ERATO 1 reads the data produced by the equilibrium code and maps them into ERATO form. It also prepares the analytic extended Solovev equilibrium [9] to be used as a test case. ERATO 2 calculates the vacuum contribution to the energy ~W. ERATO 3 builds up the matrices A(potential energy) and B(kinetic energy). The eigenvalue problem Ax = w2Bx is treated by ERATO 4 which makes uses of the subroutines of the eigenvalue package HYMNIABLOCK [10] included as part of the ERATO code. The eigenvector is analysed in ERATO 5. All 5 programs call subroutines of the CDC OLYMPUS package [16] which are also included. Subroutine maps are shown in figs. 11—16, calls to OLYMPUS subroutines being omitted. Tables 12—16 give lists of subroutines. The structure of ERATO 1 is similar to that of an OLYMPUS code but the order in which the initialization subroutines are called is not identical. The main programs ERATO 2—5 are dummies, and the entire code can readily be combined into a single program if sufficient storage is available. 2. Equilibrium 2.1. General equations In axisymmetric geometry, the magnetohydrodynamic equilibrium is given by VpjXB,
VXB=f,
17•B=0.
(2.1)
Writing the magnetic field as B=~ê~c,_±ê~,XV~,Ll,
(2.2)
we combine eqs. (2.1) and (2.2) into a second-order differential equation for /‘:
a /1 ~.,\ a2~ ~ -i--) = —j~r, where the toroidal current density4, is *
i/i
~--~—~-
dp
1 dT2
(2.3)
(2.4)
R. Gruber et al. / ERA TO stability code
326
A polar coordinate system (r, z, p) is used. The two functions p(i,Li) and T(i,ti) correspond to the plasma pressure and the flux of the poloidal plasma current, respectively. They are functions of ~tionly and can be arbitrarily specified. In order to prevent surface currents we require j~,(~i= t,115) = 0, where ~ is the flux at the plasma surface. The safety factor q (i~) is defined by dl.
(2.5)
The integration path follows a constant flux surface in the poloidal direction and B~is the poloidal magnetic field component given by 2B~= ViP2 = (~)2 + (~iP)2 (2.6) T
2.2. Solovev equilibrium A very simple polynomial solution has been found by Solovev [9] and is characterized by 4iP
— — —
2)(i~i iP~) 5(i a2E2R2 +E
~2
—
— —
~2 °
—
4iP~b2 a2R2E~i~,
27 ( ) .
leading to iP 2+b2 ~‘a2R4[5[r E2
2—R2)2 Z 2(r
.
28 ( . )
In expression (2.8) the quantity E denotes the ellipticity, R is the distance of the magnetic axis from the toroidal axis and a is the inverse aspect ratio. In this equilibrium_the toroidal current density j~,is flat. Note that on the symmetry plane z = 0 the plasma is in the range R.,~/f = 2a ~r
3. Stability problem 3.1. General axisymmetric problem ERATO solves the linearized ideal MHD equations [12] in variational form [4]: ~L=~(Wp+WV+Wk)=O.
(3.1)
W~,is the potential energy of the plasma: ~
fdr
VX(~XB)+(v~)(fX~)]2+yp(V .~)2 —2(v .~)2(jXv)(B~V)v},
(3.2)
R. Gruber et al. /ERATO stability code
327
where ~ is the displacement of the fluid elements and v the unit vector which points outwards normal to a p constant surface. W~,denotes the vacuum energy
wv=& fdr~v-
(iPViP)}.
=
(3.3)
0 is the scalar potential in the vacuum region. In the kinetic energy Wk
=
~-
f
(3.4)
drp(~-~)
p denotes the mass density of the plasma. One boundary condition demands that the component of the magnetic field perpendicular to each conductor c~vanishes =0. The continuity condition through the plasma surface (s) is written
(3.5)
Cj
[V X (~ X B)]~ (3.6) This general time and space formulation is reduced to the axisymmetric case by imposing homogeneity in the toroidal direction ~ (see fig. 1). This coordinate ~ as well as the time become ignorable variables after a Fourier analysis on the displacement vector =
~(iP~ x~p, t) = E~(iP,X)en~(~t).
(3.7)
Here n is the wave number in the toroidal direction and ~2 becomes the eigenvalue of problem (3.1). In that formulation, ~2 > 0 describes an oscillatory and stable solution, ~2 = 0 a marginal situation, and w2 <0 means that the equilibrium is unstable. z z
:___
a
=::::
-__
____
plasma surface
-2a AB
(Rest—i) a
Fig. 1. (s, x) coordinate system in plasma region (a) and (p, 0) in vacuum region (b).
328
R. Gruber et al. /ERATO stability code
3.2. Choice of independent and dependent variables Since in ideal MHD the fluid elements are frozen to the flux surfaces it is crucial to take 0 or any function of 0 as an independent variable. We found it convenient to choose as independent “radial” variable s =
where i,(i~denotes the toroidal magnetic flux at the plasma surface. The “poloidal angle” x is chosen in such a way that the Jacobian becomes J = qr2/T. In the infinite aspect ratio limit, J = J(0) only. This angle x has been used by the Princeton group [6] and proves to be convenient, since dX=q~ 11p1d1=j~-~1 ,
(3.8)
JB- Vfrr~/_+inqfEEqF(f)
has only a multiplicative equilibrium factor q = q(i,Li) in F. As dependent variables we choose the contravariant components 2,V 2siPs~x(rIViPI Y+t~3~~Xand Y=2s0 X=JIVi~’I~~/r 5~/Tr, —
parallel to the magnetic field. The quantity i3,~=
(3.9)
(axIa~)~ measures the non-orthogonality of the coordinate system.
3.3. Potentialand kinetic energy in axisymmetric geometry With the Fourier expansion (3.7) the variational form (3.1) becomes 2~K0. (3.10) 6L6Wp+~Wvw Introducing the independent variables (s, x) and the dependent variables (X, V, Y), the potential energy W~, becomes W~=1~ff~~~dx[aIIiI2 +b11 2 +c11 2 +d11 2 21 31 41
—
elI 5
(3.11)
2],
with
iax
~ ax av
11F(X)su——+inX, q8x
~
13 14
=
V ~ .
r1
2x aras
T 1 r
12=rVT—i=—+—, r as ax
a(Tx/q) + 2j~ir as lViPls ar2v q2 2
ax
~,
F(r Y),
+—~-
r
X,
I~=
(3.12)
and a
2q2i,lir4/J3IViPI2,
b
T2r2/20
2r2I2iP 5J,
c
IV’P1
4’y~/2iP 5J,
dr
8J,
4 4 ~-(lnlViPDn _.~~~(lnr)n]. (3.13) e_.±[j_ 12+ We have chosen this form in order to be able to treat the limiting case of a currentless plasma. This formulation has also analytic advantages since there is only one negative quadratic term which in a marginal situation must cancel the four positive quadratic terms. In a tokamak ordering (e = IViJ.i IT << 1), for example, the second term seems to dominate. Whenever ‘2 is not small, this term defmes the fast modes. For the other two classes of the spectrum this term must be small (12 ~ a for the slow waves and 12 ~ ~2 for an Alfvèn or a kink). Another advantage of (3.12) consists in the presence of Y only in one single, the fourth, term.
R.
Gruber eta!. / ERA TO stability code
329
The “kinetic energy” becomes
411 ~dX{fIXl2 +g(V+Y—f3~XI2+hlYI2},
(3.14)
Kr~ with
f= 2p~iTr2/qlV~il2, g = pIV~I2qr4/2T~
4Tq/2~ 5, h
pr
5.
(3.15)
In the plasma region we have to impose a regularity condition at the magnetic axis X(s0)0.
(3.16)
When we are only interested in internal modes we can put the wall right at the plasma surface. The vacuum contribution then tends to infinity and is replaced by the fixed boundary condition X(s1)0.
(3.17)
In addition we only look for symmetric displacements which are given by the condition (X, Y, V)(x) = (_X*, y* V*)(21r where
*
—
x),
(3.18)
denotes the complex conjugate.
3.4. Vacuum energy in axisymmetric geometry In the vacuum region the scalar potential 0 which appears in the vacuum energy (3.3) satisfies 20=0.
(3.19)
V By use of the vector identity =
V• (OViP)
—
0V20
(3.20)
and Gauss—Ostrogradski we obtain
(3.21)
Wv~ifO~ds+~~ fords,
where p and c represent the surface of the plasma and of the external conductors, respectively. For n = 0 and for a symmetric displacement, the compression of the field lines in the inner region of the torus leads to additional terms [131. On each conductor c• the boundary condition implies = 0. The continuity condition on the plasma surface (s) is
[VX(~XB)] Then W~becomes W~ —T ~
f
v=~ ~
r0~F(~ d~.
F(~.
(3.22)
(3.23) (3.24)
R. Gruber et a!.
330
/ ERA TO stability code
Including the boundary condition (3.5), we can use Green’s theorem 0(x)
=
~L f[o(x’ 2ir S
aG(x,x’) a~
G(x, x’)
—
~~11 ±~ ~ a~ ds’
—
fO(x’) aG(x,x’) ds’
2ir
(3.25)
ci
with G(x,x’) = 1/lx x’l, where x’ and x are the integration variable and a point, both on the boundary surfaces, respectively. In the numerical treatment it will be convenient to use polar coordinates (p, 0, p) to integrate (3.24). In this system: dz 2 i=~i~_, ds=lrd~,d0=rdl=JlViPldx, 3 l/dz a dr 3\ i/7dr’\2 —
~
+(_),
r= R +p cos0,
Ix —x’I
p sin0,
z
=
[(r
—
r’)2
+ (z
—
z’)2 + 4rr’ sin2 (p12)]112.
(3.26)
After having Fourier analyzed 0 and 3çb/&v on the surfaces 0(0, ‘p)l
5~= 0(O)I~,~e’~, 30(O,~ =N(0)e~’~’
(3.27) (3.28)
ir av and defined the p-integrals G~(O,0’) e1’~
2ir
,
G~n’(O,0’)=f ~
dp,
(3.29)
IX~X~’~
we obtain the two following coupled equations: 0~(0)=
~—fdO’ r’v’
0~(0)=±
f
do’ r’v’
V’G~’0~(0’)
V’G~ 5’Ø5(0’)
—-~—f dO’ r’v’- V’G~’0~(0’)
(3.30)
—~fdO’ r’v’
(3.31)
do’ G~’N(0’)
—
1f
dO’ G~’N(O’)
V’G~c’0c(O’)
—
with v
•
V
,~ dz’ 1aG~B’\
GAB’
—
~
\
a,~’J~,
—
dr’faG~B’~ az’ ~
ã~k
(3.32)
and A and B can both either be on the plasma surface (s) or on the conductor (c). The prime is set on the integration variable. In the present version of the code we restrict ourselves to one conducting wall only. 3.5. Additional n
=
0 terms
In the previous treatment of the vacuum contribution, we assumed that for the perturbed magnetic field in the vacuum region: jIiB. dl=~ViP- dl=0.
(3.33)
Condition (3.33) is always true for displacements with n 0, and for n = 0 and no m = 0 contribution. However, for n = 0 and m = 0, there are other terms which have to be added to W.,. due to a compression of the average
R. Gruber et al. / ERA TO stability code
331
magnetic field by the displacement. The perturbed magnetic field liB in the vacuum region can be written
6B=zi~+z2X~
(3.34)
,
where Z1 and Z2 are constants to be calculated, a fulfils A*ct=0
ac =0
a51,
(3.35)
and 0 denotes the scalar potential which has been treated before. The vacuum energy then becomes 2 ÷z~ (e~~X Va)2 W~= 4 J’dr
[z~(~)
(3.36)
Note that all the coupling terms disappear due to (3.33). In this section we calculate the first two terms of (3.36) for which a closed line integral in the vacuum does not vanish in all the cases of an axisymmetric mode. To compute the constants Z 1 and Z2 we form the integrals
f (dsp XA)
Ii
12 =
•~,
~
f(ds~XA)
(3.37)
where A is the vector potential corresponding to SB defined in (3.34). On the plasma surface d.c = rdp dI X ê~,. Using the boundary condition
(3.38)
(vXA)~=—~,•B, the integrals (3.37) become 11
=j~(dlXv)
B,
‘2
(3.39)
=j~,p(dlXv). (BX Va).
The integrals can also be converted into volume integrals over the vacuum region Ii=fdrV.(AX~!~L)’*’Zifdr3~,
I2fdrV.(AX~~)Z2J’dri-!,~_.
(3.40)
Here we have used the relation (3.33). By equating the expressions (3.39) and (3.40) we obtain the constantsZ1 and z2 in terms of ~ on the surface. By substitution into the expression (3.36), the final vacuum energy becomes 2/fdr-~ + [j~@1xv). (BX Va)]/ fdr~!~_}. (3.41) W~= 4([1~(d!x v) B]
Knowing that —4n~’Gk~’ is the Green’s function related to the operator
1=~
~-~-~-
_~-_J
v’. V’(—4rr’G~~’) —
~-
v’~V’al’(—4n-’G~ 5’)+
we can write Green’s theorem
~
~J— I
V’al’(—4rr’G~’)
(3.42)
~—v’. V’al’(—4,r’G~’)
(3.43)
‘•
on the plasma surface and 0
=
i-_f
~-
v’- V’(—4,r’G~~’)
on the conducting wall.
—
f
~
v’. V’~’(—4n~’G~’) ~
f
332
R. Gruber et al.
/ ERA TO stability code
In (3.30) and (3.31) we knew ViP and had to calculate
0. Here, we know a and have to calculate Va.
4. Numerical approach 4.1. The marginal point of the continuous spectrum In the 1D stability code THALIA [14] we had to choose our basis functions such that we were able to reach the marginal situation (w2 = 0). This spectral point was characterized by the condition V•~=0. (4.1) Spectral pollution [5] then disappeared and the stable part of the spectrum proved to be well separated from the unstable one. We adopt the same idea to solve the 2D axisymmetric stability problem and start to calculate the conditions defining the marginal situation in toroidal geometry. They are given by the three identities [15]:
(4.2)
F(V)=-~-+inV’0,
q ax
(4.3)
2
r
3s
r
—
3x
+~j- F(Yr2)~’0.
(4.4)
r
Eq. (4.2) suggests a Fourier analysis in x~[6], eq. (4.3) imposes the choice of basis functions in s and eq. (4.4) simply defines Y. Another possibility to fulfil eqs. (4.2)—(4.4) is the use of finite hybrid elements [7]. 4.2. The finite hybrid elements These involve extending the variational formulation and making a finite element approach in s and of solving ~L(X, V, Y)=0,
x~Instead (4.5)
we vary over 7 variables and solve 6L(~-~-_--, ~(2)~~T~
~
j~(2)
~
y(2)) =
0,
(4.6)
restricted by the evident identities ~1)
=
=
~
ji(l) =
v~2~ y(l)
=
y(2)
(4.7)
The original problem is reobtained when the identities, eqs. (4.7), are substituted back into eq. (4.6). In our approach, however, we vary the Lagrangian over all 7 variables ~ ~ ~(3) ji(1) ~(2) Y(’) and y(2), and do not impose the identities in eqs. (4.7) everywhere, but only on specific points in each mesh cell. To prepare the discretisation of the problem, we rewrite the constraints (4.7) in the equivalent forms urn
—~-
~
f(x(2)
f(v(2)
—
X~’~) dr = 0,
V~E ~l,
—
v0)) dr = 0,
V~E ~
lirn
-~--
~
f(x12)
f(y2
—
x~3~)dr = 0,
Vi~E ~2,
—
Y(1)) dr = 0,
V~E ~l,
(4.8)
R. Gruber et al.
/ ERA TO stability code
333
where ~l denotes the plasma domain. The conditions (4.2) to (4.4) then become:
av~1~
(4.9) =
+
0, 3lnr2
~)
~2)
(4.10) + ~2)(inq
(4.11)
+3~2)+~~~=o.
4.3. Choice of the mesh In order to discretise the problem we cover the domain ~l with a rectangular mesh. In the radial direction we choose N~intervals (the radial mesh point counter i varies from 1 to N 1 + 1) and in the poloidal direction we choose I / ~ N,, + 1. All the mesh cells can have an arbitrary size. In order to be able to impose the symmetry condition (3.18) we add two half mesh cells in the plane ir ~ x ~ 2ir (see fig. 2). Then we can confine the solution to the half plane 0 ~ x ~ ir only. We now satisfy the constraints (4.8) by identifying ~ as the size of a mesh cell L~,j,such that 2)—X~3~)dsdx0,Vi,/, I (X~ —~‘~)dsd~0, Vi,j, (X( ii ~ij ii ‘~
±
f f
~—
(V(~~ —
V’~)
ds dx = 0,
Vi, j,
_L
ii ~ij
f f (y(2)
—
Y(1)) d.c
dx = 0,
Vi, j.
(4.12)
‘I
These conditions do not uniquely determine the elements, so that we still can choose the order. We make the simplest choice possible by requiring that each argument in eq. (4.6) be constant on each mesh cell. Calling d 1 the basis function1’for a11 the basisthefunction for can ~ beg1+1~2,~ the and~hi+l/2J, thebasis basisfunction functionfor for~v<2) f~ andthe y(2), expansion written: basis function for V~’~ and Y~ NS+1 Nx+l
NS+1 Nx+l = ~
~
2~ = I~
X 11d11,
i=1
j=1
x~
i=1
NSfl Nx+l
~
j=1
2)
1=1 j1
(V) g~+112,1, Y 1+1/2,1
11j~1,
x~
1~11,
NSN+1
(V)(l) =
x
3~ =
x,
~
~
= N~N+1
(V)( i=i 1=1
(1
i=1
Y
1+1/2,1
j=1
h
1+1~2,1,
(4.13)
where the positions of the nodal values in a mesh ate all shown in figs 3 and the shapes of the corresponding basis functions are shown in fig. 4. It is easily seen that the expansion (4.13) together with these basis functions fulfils 2) reduce to the requirements (4.12). With this expansion the 7 unknowns ~ ~ ~(3) 1/(1) 1/(2) y(l) and y( the 3 unknowns per mesh cell X. 1, V1÷112,1and Y,~112,1. Note that d11(x) and ~ are functions linear in x and piecewise constant ins. The x derivatives of such elements are piecewise constant. A similar roof function linear in s and constant in x is used for f~1(s)such that its s-derivative is piecewise constant. The bases for f~1and h1+1~2,~ are piecewise constant. The basis functions for the radial component X are non-zero over 4 mesh cells, whereas those for V and Y are non-zero over 2 mesh cells only. The conditions (4.9)—(4.1 I) then can be fulfilled when they have to be fulfilled. We can now easily write each term in the Lagrangian (4.6) at the mesh centre:
R. Gruber eta!.
334
I ERA TO stability code
Ix
x
j 6
— .~_.
x~,
——
1,J.1
VL,lh 4
—
3
—
j.l
Nx5 N5=4
2
s
0—1
2
3
45
VL1/2j
Fig. 3. Position of the nodal values in a mesh cell for the “finite hybrid elements”. The physical values are situated in the centre of the mesh cell (o).
-
Fig. 2. Discretization of the domain l~.
ax~+112,~+1,2 x~+1,1+~ +x1,1+1 — ~ 2(Xj+i —xi)
X~,i.1
—
ax
x~2~
_x~~÷x
1,1+1 +X~+1,1+X1+1,1+5 4
—
_____________
—
v
v
Ai+l,j+l
Af~1,/— A1,j+l
2(s,+i av~+~/2,J+l,2 = V~+112,1+1— —
__________=
Yi+l/2,f+l
—
—
V1+l/2,/ x~
1/(2)
=
1+1/2,/+1/2
‘
2)
y(
~‘i+l/2,/
—
Au
s~)
=
V1+j/2,/+l 2 Y
+ V1+112,J
1+i/2,/ +
Yi+l/2,i+l
(4 14)
i+1/2,j+1/2 2 These relations can be considered as centred finite differences applied to a variational problem. Since in the Lagrangian, eq. (4.6), all the functions and their derivatives are piecewise constant, the equilibrium is also taken to be piecewise constant. This is done by replacing all coefficients coming from the equilibrium by their values at the centre of each cell. With (4.14), the boundary conditions, (3.16), (3.17), take the simple form —
x~,I=0, xJ%J1f1,J=o,
‘
l~/~N~+1,
(4.15)
and the symmetry conditions (3.18) in x yield X~’—X~, X~1=X~2, 1~i~N1+l,
~ = ~ ~ 1 = ~+l/2,2 V~÷1121 =—V~+112,2, Y~+l/2,l= Y~+1~,2
1
~ i ~N5.
(4.16)
We can now verify that all the three conditions (4.9)—(4.1 1), being piecewise constant, can indeed vanish over each mesh cell.
R. Gruber et al.
/ ERA TO stability code
7
j+1
~
335
d
e
j
.1~
~
1 ~
x
j+1 ~
~
+
g
~
1/2
Fig. 4. Shape of the basis functions d 2~, of f~ 3), ofg,~ 11(x)0)for and they(l) expansion and ofh~+jI of XW, of e11 for the expansion of 2)ax( nd y(2)~ 1for the expansion of x( 1j2,fom the expansion of V 2,jfor the expansion of v~
336
R. Gruber eta!. / ERA TO stability code
Table 1 The 5 ERATO main programs Program
Contents
ERATO 1
initializes, chooses equilibrium, interface equilibrium, stability, Mercier’s criterion, ballooning criterion
ERATO 2
calculates vacuum energy contribution
ERATO 3
calculates potential (A) and kinetic energy (B) matrices
ERATO 4 ERATO 5
solves eigenvalue problem Ax = o.~’ printing, plotting and display of the eigenvalue, the eigenfunction and other diagnostics
2Bx
4.4. Spectral shift and its elimination Compared with the usual finite element method [8] the additional freedom introduced by the enlargement of the functional space leads to a lower energy level. This helps to separate better the weakly growing unstable modes from the stable continuous region. However, for the Solovev equilibrium we observe a spectral shift [15] due to the numerical approximation of eq. (4.9). This spectral shift can mostly be absorbed by a small modification of nq. For numerical purposes it is convenient to impose the shift correction on the wave number n. This means that instead of calculating with an integer n value we introduce such that ~‘
(4.17)
i~n/{l —a(nq/N~)2—I3(nq/N~)4]. For a Solovev equilibrium we found a = 0.822 and 13 become smoother.
=
0.142. The convergence properties for fixed nq then
Table 2 Disk units Unit
Disk file name
Created
4(MEQ) 5 (NIN) 6(NPRINT) 7(NDA)
TAPE 4 INPUT OUTPUT TA
ERATO ERATO ERATO ERATO
8(NSAVE)
1 1 1—5 3
Used
Contents
ERATO 3,5 ERATO 1 ERATO 1—5 ERATO 4
mapped equilibrium read NAMELIST from card line printer unit matrix A NAMELIST on disk matrix B matrix LDLT numerical equilibrium scratch unit shifted A eigenfunction plotting quantities scratch unit surface equilibrium vacuum matrix
OUT
ERATO 1
ERATO 2—S
9(NDB) 10(NDLT) 11 (NDORY) 14~NSCRTC) 14(NDS)
TB TAPE 10 DORY TAPE 14 TAPE 14
16(NDES) 1 7(NDMAP) 17(NVAC)
TAPE 16 TAPE 17
ERATO 3 ERATO 4 equilibrium ERATO 3 ERATO 4 ERATO 4 ERATO 1 ERATO 1 ERATO 1 ERATO 2
ERATO 4 ERATO 4 ERATO 1 ERATO 3 ERATO 4 ERATO 5 ERATO 5 ERATO 1 ERATO 2 ERATO 3
R. Gruber et a!. /ERATO stability code Table 3 Common blocks No.
Name
Used by
Contents
C2.1 C2.2 C2.3 C2.4 C2.6 C3.1 C3 .2 C3.3 C4.l C5.1 C9.1 C9.2 C9.3 C9.5 N
COMPHY COMEQU COMEIG COMAUX COMVEC COMESH
ERATO 1—5 ERATO 1,3,5 ERATO 1,4 ERATO 1,3 ERATO 4 ERATO 1,2,3,5 ERATO 1—S ERATO 4 ERATO 1—5 ERATO 1—5 ERATO 1 ERATO 2 ERATO 3,4 ERATO 5 ERATO 1—5
physical quantities equilibrium eigenvalue constants auxiliary quantities eigenvectors mesh numerical quantities eigenvalue constants control values disk units mapping vacuum matrices A and B plotting quantities namelist
COMNUM COMIVI COMCON COMOUT COMMAP COMVID COMMTR COMPLO NEWRUN
5. Organization and normalization The ERATO code consists of five programs ERATO 1—5 which are explained in table 1. In addition, three external packages, namely: the numerical equilibrium (actually we use the ORNL code [11]), giving all the equilibrium and geometrical quantities to the stability calculations, some utility routines from the OLYMPUS library [16], the eigenvalue solver HYMNIABLOCK [10], are employed. All these codes are linked together by means of disk files which are described in table 2. All the disk input—output is performed in the subroutines IODSK 1—5. Table 3 contains a list of all the COMMON blocks used in OLYMPUS, HYMNIABLOCK and ERATO 1—5. —
—
—
Table 4 Normalization Position
rph = RTE
Displacement
tph = RtE
Fields
Bph = B0BE
Mass density Pressure Frequency
Pph = POPE 2o/12o)PE ‘~ph= (B ~ = (B~/(poR2~o))w~
Current
‘ph = (BORI,2o)JE
Notes: R is the real distance of the magnetic axis from the symmetry axis, B 0 denotes the toroidal magnetic field at the magnetic axis, P0 is the mass density at the magnetic axis, ‘E is the total current which is calculated and printed out in ERATO I (see section 6.9), eqs. (6.21).
337
338
R. Gruberetal. /ERATO stability code
At the end of section 12 we give a subroutine map of each main program ERATO 1—5 and of the complete code. Table 4 contains the normalization of the equilibrium quantities, the eigenvalue w2 and the eigenvector ~. Table 5 Namelist NEWRUN Quantity
Type
Default value
ALARG ALO ANGLE(16) ARROW ASPCT B2R2 CPSRF CST ELLIPT EPSCON EPSMAC GAMMA QIAXE
R R RA R R R R R R R R R R
0.0 —1.5 0.0 1.0 1/3 0
REXT SCALE WIDTH WNL WNTORE ITEST LENGTH MEQ NAN NCHI NCOLMN NDA NDB NDES NDLT NDS NFIG NIN NITMAX NMESH NPRNT NPSI NR NSAVE NSCRTC NSHIFT NSUR NTCASE NTURN NVAC NZ NLTORE NLY
R R R R R I I 1 I I I I I I I I 1 I I I I 1 I I I I I I I I 1 L L
4 10~ i0’2 5/3 3.3333 0.0 1.0 1.0E5 —2.0 1 4(Nx 4 ~JN5 15 8(Nx 7 9 16 10 14 0 5 10 0 6 14 65 8 14 0 15 2 10 17 33 T T
+
l)(8Nx
+
1)
+
9)
paper size for plotting eigenvalue guess angles at which displacement plots are performed length of arrow inverse aspect ratio a b2R2 for analytic equilibrium set in ERATO 1 set in ERATO 1 ellipticity squared E2 accuracy of eigenvector machine accuracy specific heat ratio inverse q-value on axis 1/q 0 defines conducting wall scaling factor parameter for non-equidistant x mesh dummy parameter toroidal wave number n CDC plotting parameter defined in AUXVAL channel number for equilibrium defines region to fit on number of poloidal intervalsNx defined in AUXVAL channel number for A channel number for B channel number for plotting channel number for decomposed A scratch unit number of displacement plots input channel number maximum number of iterations choice of mesh output channel number number of radial intervals N~ number of r intervals in equilibrium channel number for NEWRUN scratch unit shift in q plasma surface for plot selects the case number of turns for ballooning criterion channel number for vacuum number of Z intervals in equilibrium T = toroidal geometry T = compressible case,F= incompressible case
R. Gruber et al.
/ ERA TO stability code
339
6. ERATO I
6.1. Initialization ERATO 1 first reads the 4 label input cards. In CLEAR all the COMMON variables should be put to zero. PRESET then sets the first test case later discussed in section 11. The NAMELIST input (see table 5) is read in DATA and in AUXVAL we calculate some matrix sizes used in later main programs. Then FORMEQ selects all the cases defined by the NAMELIST quantity NTCASE, described in table 6. For NTCASE = 3 or 4 numerical equilibrium quantities must be transferred via disk. This is foreseen in IODSK1 fork = 6 (table 7) which reads the number of mesh points in rand z, the equidistant rand z mesh 0(r, z), and, depending on the numerical equilibrium code, constants to defme the arbitrary functions p(~) and T(0). For our equilibrium code this statement is
READ(NDORY) NR,NZ,((CPSEQ(JR,JZ),JR 1 ,NR),JZ=1 ,~NZ),(EQR(JR),JR1,NR), (EQZ(JZ),JZ= 1 ,NZ),CPSRF ,(APP(J),J=1,1 0),(ATTP(J),J=1 ,10)
(6.1)
The way one reads NR, NZ, CPSEQ, EQR, EQZ, CPSRF, APP and ATTP is free; it is only necessary to be consistent with what has been written on disk in the equilibrium code. The two arbitrary functions p(~) and T(0) have to be defined in TANDP and must agree with the functions chosen in the equilibrium code. They have also to be normalised such that they correspond to eq. (2.3). In our equilibrium code the chosen functions are
TT’(O) =
=~
~
ATTP(J) ~‘.
(6.2)
The quantity Po is defined by the condition p(~ = 0~) = 0. The normalization with respect to R and T(0 is performed in subprogram NOREPT. The plasma surface is defined by passing the value i~ion the surface, i.e. CPSRF(tL’5) through (6.1).
0)
6.2. Choice of the mesh The existence of singular nq = m (m integer) surfaces which are concentrated at the outer plasma region requires a careful choice of the mesh in the s-direction. This is done whenever NMESH * 0. First we go through the whole mapping procedure with an equidistant s-mesh in order to count the number NSING of singular surfaces for which nq = m.
Table 7 IODSKI(k) k Table 6 NTCASE ______
21
rewinds all disk files writes NEWRUN on NSAVE
NTCASE ______________________________________________________
3 4
writes EQ on MEQ writes plot quantities on NDES
1 2
analytic Solovev equilibrium analytic equilibrium through interface
3
numerical equilibrium through interface fitting of the plasma surface
5 6 7
writes vacuum quantities on NVAC reads equilibrium from NDORY not used writes surface and labels for plot
4
____________________
8
340
R. Gruber eta!.
/ ERA TO stability code
In a second step, the distribution of the mesh points in the s-direction is controlled through a monotonic weight function W(s) satisfying W(0) = 0 and W(1) = 1. The locations of the mesh points s~,i = 1, N5 + 1, are ...,
given by s~= W’((i
—
1)/N5), ‘c/j,
(6.3)
1 (x) is the inverse weight function.
where s= W x Designating
1, / = 1, NSING to be the locations of the NSING singular surfaces lying within the plasma, the density of mesh points is chosen to be 3-i-NSING NSING 1 2+~ (6.4) ds3(1+NSING)~ 1=1 (s—x1) For the width L~we use = [ndq/ds];~.The normalization constant a is determined by W(l) = 1. The weight function then becomes
dW
3+NSING W(s) = 3(1 + NSING)5
NSING +
a
j=1
arctg
s—x~ +
arctg
x ~—
,
(6.5)
where 2NSING a=
1NSING ____________
3(NSING+1)i
‘
1=1
arct g
1
—~.
+
arct _.1~ g~
66
When NSING = 0 the mesh remains equidistant. The constant term in the density (6.4) prevents too strong a concentration of points around the singular surfaces. The formula chosen does not pretend to be the best one. We have sometimes multiplied the ~j’s with a factor ~ to put more points around the singular surfaces. In the angular direction x we have chosen an equidistant mesh. It has appeared that it is sometimes advantageous to concentrate the x’~near the outer plasma region, i.e. around x = 0. This is specifically true when a ballooning type of mode is to be represented. The user is free to modify the subroutine MESH. 6.3. Analytic treatment
of
the axis
2z2) around the magnetic axis we perform a 10 parameter fit To obtain an adequate knowledge of ~/i(r ~L’ —~f~ +f 2 +f 4 +f 6 +f 2 +f 4 +fjoz6 +f 2z2 +f 4z2 +f 2z4 (6.7) 3r 5r 7r 8z 9z 12r 14r 16r over (NANR—NANL)*NANZ equilibrium points. This region around the magnetic axis is defmed by NAN which denotes the number of s surfaces which should lie in the “analytic region”. The default value is NAN = ~ The magnetic axis R = RAXIS is found by solving ,
aiP/ar(z=0)=0=f
2+3f 4. 7r
(6.8)
3+2f5r 6.4. Mapping
We assume that the numerical code furnishes ~ (r, z) where the r and the z mesh are equidistant. For the stabi-
lity code we have to invert the information, i.e. for a given tji we have to find the corresponding r and z coordi-
nates. This is performed in subroutine INTFAC. We start our treatment at the z = 0 axis in the outer region of the plasma and first fmd r(iP, x = 0) by a linear interpolation between the two neighbouring radial mesh points. By means of the subprogram ADVANC we
R. Gruber etal.
/ ERA TO stability code
341
~
r(s,X=lr)
R
r(s.X~o)
Fig. 5. The 6 points which define the polynomial (6.9) in the mesh cell (r
1
advance in the positive x direction giving for each mesh cell the two ingoing points (nos. 2 and 4 in fig. 5) between which ~,1ilies, and searching the two outgoing points (nos. 4 and 5 in fig. 5). The spatial of the 1 = 2 anddependence z2 constant line is, for each mesh cell, approximated by a third order polynomial in r ~i(r2, z2)”a 2 +a 4 +a 2 +a 4 +a 2z2. (6.9) 1 +a2r 3r 4z 5z 6r To determine these 6 parameters a 1to a6, this polynomial goes through all the 4 edge points (nos. 2 to 5 in fig. 5) of the mesh cell. The two additional points are chosen to lie inside the surface and so that there are at least 3 distinct rand z coordinates. When we cross the z = 0 plane again, the angle x (3.8) must be equal to ir which defines the q-value of this =~/‘ surface (2.5). The ~,1iderivative of q is given by the periodicity condition on the non-ortho13~(j3,~(O) I3~(7r) = 0) gonality 1 .-i’r i /~ç~~2 dj3~= dx ~j-~- i~i~ (6.10) -~~-+
1.1 ~-
—
~—~-~----
— —
which leads to
dqq t’d
j
/~,r
X
1 dT +
1
Iv~I2~
611 .
(.
This simple procedure is very accurate as long as we are not too close to the magnetic axis which is achieved by the 10 parameter fit (6.7). 6.5. The normalized equilibrium quantities In the subprograms EQANAL (analytic treatment) and in NUMEQU (numerical mapping), we fill up the array EQ(20,NCHI) which contains all information about the equilibrium for an s constant surface. The 20 quantities in EQ per mesh cell are explained in table 8. All the EQ(20,Nx) arrays are stored on disk file MEQ = 4. There are N 5 such arrays on that file. In all the quantities defined in table 8 we assume that the flux function Ton the axis is normalized to 1 and that the distance the toroidal axismeasured r is normalized 2 thatfrom is obtained is then in unitstoofthe distance ofthe magnetic axis R (see table 4). The eigenvalue w w~B2 2, (6.12) 1]popR where BT is the toroidal magnetic field at r R.
342
R. Gruber eta!.
/ ERA TO stability code
Table 8
EQ(k,i),j1,Nx
-_______________________________________________________
k 1
s-coordinate of lower left hand edge
2
x-coordinate of lower left hand edge
3
s-coordinate of upper right hand edge
4
5,6
k-coordinate of upper right hand edge (s, x) coordinate
7
mass density p
8
plasma pressure ‘yp/qo~l
9
flux
10
free
11
q/qo
12 13
poloidal magnetic field ~pr non-orthogonality i3~= (axIas)~
14
distance from toroidal axis r2
15
(a lnr2/as)~ (a In r2/ax)
16
5
function T
2fqolV~ 2
5
2
17
~
18
K=—I
slVVd
T ds \q q
2~r ~2
/
dp a
a
___~_~+_~!.
—(In IV~I)~———(Inr)
qoLlV~PI r a~p
dXL’
2/J=T/q 19 20
a~
0
rfree
6.6. Scaling law By means of a scaling factor a (in the NAMELIST input SCALES), it is possible to transform one equilibrium into another one [23]. Calculating one numerical equilibrium makes it possible to perform a whole stability study in q 0 as a parameter. At the same time we set to 1 the toroidal magnetic field at the magnetic axis T0 = 1. The scaling law is given by 2p, ~—~-al/i, ~~-+aB~, ~-~‘a aT~, D2 T~+ a(T2 Tg), ~ -~q2~ (6.13)
4-+a/~,
~
-
-~
6. 7. Mercier’s criterion
Mercier’s criterion [26] can be written in the form M=~((
1)2
2(l+q)
_P’T~(Ci+C l+q
2)÷P1Cs(C4+T2Ci)+p12{~(~_C1C3)_C4(C1+2C2÷C3)]}~0. (6.14)
R. Gruberetal. /ERATO stability code
For a given ~-surface, the integrals C1, C2, C3, C4 and C5, with 2/T, J’qr
343
(6.15)
become
q
~
d~
Ci~J
—
C
~
q
~
(i2
—
R2)2 d~
2~J
2
22
(r —R) dx
C42ir~,
C5
rf
(6.16)
~
This form has been incorporated in subroutine MERCIR which also uses INTGRA for the ~-integration and QUADRA for the parabolic fit. As a result we getM~0 for a stable situation, while M <0 predicts an unstable localized mode. 6. & Ballooning criterion On each constants surface we test in subroutine BALOON if the modified potential energy [17,21] = ~-fdx-~-
(A
I~f2~+BXoI+ CIX0I2J
(6.17)
is positive definite (ballooning stable) or not (ballooning unstable). The quantitiesA, B and C are given by [see (3.13) for the definition of a, b, c and e] q’(x XoY~ q / abc ( 2jVir T’\2 C=(b+c4~\sIVlpI2+~~:;/ —e.
A
=
a
+ —~--I’13 +
—
B=
‘
bc
(13
(b+c)A\”~
+ ~‘(x
—
q
T2/’ Xo)~(2/iPr + T’\ /~~V~p(
(6.18)
These expressions have been obtained after expanding X and V for large n as
~
(x 0
v= (v0 ~
~
(6.19)
As periodicity condition we impose
Xo(~= —lir)
=
(6.20)
Xo(~= iir) = 0,
choosing ito be a large number (up to 100). The printed array BALLCRIT contains the number of negative eigenvalues of the 5 W~,matrix where the values 0 mean stability and the values >0 instability. Close to the magnetic axis 1 and inprecise information can be produced. However, we know that around the axis MERCIER and BALLCRIT are equal. -~
6.9. Definition of the betas Let us define the integrals (in BETAST) 2, I~ di~d~p, I~2=JfJd~d~p
~ff~
1B2
fJ”1 d~d~,
344
R. Gruber et a!.
zv-]fJdiidx,
/ ERA TO stability code
IE=f~f~id;l1dx/~,,lRff~d~l1dxp.
The usual beta value (/3), the poloidal (j3~)and the Princeton beta (/3*) are given by (3 = 21p/1B2, i~1= 81rIR/I~, (3* = 2~/’l~i~/IB2.
(6.21) (6.22)
6.10. Surface fitting The calculation of a numerical equilibrium takes a lot of computing time whenever discrete conductors are prescribed. In this case, the resolution of the equilibrium is not good enough to give accurate information about stability. To eliminate this problem, we pass twice through equilibrium codes. We follow the steps (I) Calculate the numerical equilibrium prescribing positions of the wires and the currents flowing in them (ORNL code). (II) Fit the position of the plasma surface (NTCASE = 4). (III) Solve the numerical equilibrium prescribing the plasma surface as found in (II). The resolution of such a code is much better and, in most cases, enough for stability studies. In this section we described the part of ERATO 1 which performs step (II). In this step we keep the polynomial factors of eq. (6.9) for all (r, z) mesh cells through which the surface curve s = 1 goes. For practical purposes we propose to use an efficient free boundary equilibrium code including a “Buneman Solver”. Such codes enable us to resolve the equilibrium problem with sufficient precision (up to an (Nr, N~) mesh of 512 X 256) [24]. An ideal way to resolve the fixed boundary problem would be a finite element approach in which the nodal points are iteratively put in the same ~(i = constant surfaces that are required by ERATO. Such a code already exists (SELENE [25]) and will shortly be combined with ERATO. 7.ERATO2 7.1. Treatment of the singularities In this main program we calculate the vacuum energy contribution expressed in terms of the normal displacement X on the plasma surface. In eq. (3.21) we replace 0 and aO/ap by X using eqs. (3.23), (3.28), (3.30) and
(3.31).
Defining the modified elliptic integral K 0 to be
1/2 [cos~~~n~]1/2
K0~~) =
(7.1)
and
2 + (z~’ zA)2 ]/[(r’a’ + rA)2 + (Z’~’ ZA)2], [(r’B’ rA) the real and the imaginary part of the integral (3.29) become
(7.2)
G~~’ =4K
0/[(r’~’+rA)2 + (z~’ ZA)]. 0(~)(—1) The integral K 0 as well as its derivative dK0/dr! exhibit singularities, since for small n
(7.3)
Kn(~)=(_l)0~1(~)+...,
(7.4)
—
—
=
—
= (_1)0+1
—
±~+ ....
(7.5)
R. Gruber et al. / ERA TO stability code
345
The logarithmic singularities (7.4) in the source terms of the two coupled equations (3.30), (3.31) will later be treated analytically. It turns out, however, that the poles in the integrals containing u’V’G~~’ must be studied carefully and, numerically, be treated in a consistent way. Only when this is done properly are we able to move the wall continuously towards the plasma surface [22]. In our case, we treat these singularities by adding and subtracting the same types of integrable singularities, namely d0’r’v’~V’G~’= 0,
2ir
~ ~-
~—fdO’r’v’
dO’r’v’’ V’G~°5’ dO’ r’v’ V’G~~’ = —2.
f
V’G~°~’ = —1,
(7.6)
These equations mean that the solid angle of a closed surface is 0, —2ir or —4ir if the observation point is outside, on or inside this surface, respectively. Eqs. (3.30), (3.31) can then be written
~—fdO’ r’[v’. V’G’1’O(O’) v’~ ~—fdO’ G~’N(0’)_~_fdO’ r’[v’ V’G~’Oc(O’) —
2O~(O) 20~(O)= —
—
—
0~-
f
v’G~~’o~(o)],
(7.7)
dO’r’[v’. t7’G~’0~(O’) V’. V’G~° —
5’çb~(0)]
f
— ~—
V’~
dO’ G~’N(O’)
—
~—fdO’r’
[ii’.
V’ G~~’ 0~(0’)
—
V’~ V’ G~c’0c(O)].
(7.8)
All the integrals containing v’ V’G~B’are now regular. The two source terms containing G~~’ still exhibit a logarithmic singularity which we treat analytically. Instead of solving
f
5AB’ =
do’ G~~’N(O’),
(7.9)
we solve = 5reg + Sanai,
=
~
f
Sreg =
f
do’ ln[~2 + (0
—
do’(G~n’N(o’)+
O’)2],
~2 =
(p
ln[~
+
(0
—
01)2 ]N(o)J~
—p’)2/pp’,
where Sreg is a regular integral and Sanai is the part which can contain a logarithmic singularity.
(7.10)
(7.11)
346
/ ERA TO stability
R. Gruber et a!.
code
7.2. Discretisation In order to solve the system of coupled equations (7.7), (7.8) we expandN(0), Ø~(0)and O~(°) in terms of piecewise constant finite elements fI+i/2 Nx
N(O)
—
Nx
Nx
0~~O)= i0 ~ b
i=0 ~ a1+1/2 f1+i/2~
Ji+1/2~
1=0 ~
—
c f+l/2fJj+1/2~
(7.12)
integrate over 0 and end up with a set of coupled matrix equations 2b—2cAb—B•a—C•c,
0=D~b—Ea—F’c.
(7.13)
The matrix elements are
1
01+1
A11’= ~
1
~
2ITA0 1+112
dO
~
(f
dO fi+l/2
f
0j+1
ln[~2 +
2+
do’ fj+l/2 (G~s’+~Lln[~
0~
01+1
O’)2])
U
dO fi+l/2
(
~0i+1 dO fi+i/2
Ojf~
dO’ r v
V G~’
2ir
dO’f/+l/2 r’v’
I~0j
0
f 0
01+1
1
f
1
fl+i/2~1/
0~
61
1 2ir~01÷112 Ui
dO’ r’v’
2ir
dO’fj÷q 12r’v’ V’G~c’
=
—
I 0~j4~
dO f,+l/2
_____
CS
(0
J”
(0 _O’)2] dO1)~
6i—1
1
sc
fi÷1/28i/
—
01+1
f
B”= 2zM1~112 1 ~fi+l/2
dO’ fj÷112rv V’G~’
fi÷i/2
Oi+1
2sr
6/÷i
f
V’G~’
dO r v
fi+1/25i/ 0
—
2 + (0
do’f
—
I
I
~I~-,Q ~ ~,~S’j
~
01)2])
1÷112~G~5’+1 ln[~
~cs 0~
ln[~2+(O _0I)2] 01+1
Fe”
= 21r,~.Oj+1/2r 1
0/+1
dO fi+l/2 6
cc
do’).
I foi
27T
dO*fj+i/
f
2rIvI. V’G~~’fj+1/2~jf0 —
~.
r’v’•
V’G~~’ )~.
(7.14)
The expression v’’ V’G~B’in polar coordinates becomes
rag
aK01
~ __________
=
(a2~”~ + 4r’r)2
1 1a
2
(
2 p’cosO ,
a
ag
,
az’
op +~I
sin o
)
+ (p
2
—
p’)
+
2p[p
—
p
ap sin(O
—~
—
0)]}~
(7.15)
R. Gruber et al.
~ r, r’) = 2
a
=
=
/ ERA TO stability code
p)2 + 4pp’ sin2
—
~I
347
Note, that a can go to zero only when A and B both denote the plasma surface or both denote the wall. 7.3. Vacuum energy matrix Eliminating c from eqs. (7.13) we obtain b
a,
=
0
=
[A —21 —(C
—
21)F’D]’ [B—(C
(7.16)
21)F1 E}.
—
Whenever the wall goes close to the plasma (hEll << MAll) C~A—E,
D~A+c,
F~A,
E~B,
(7.17)
and Q becomes
(7.18)
Qn~E’B.
Note, that for a wall tight on the plasma all subdeterminants of Q are This implies that q.~0 which corresponds to the condition (3.17). The way that we calculate the matrices enables us to pick up this singular behaviour of Q and therefore it is possible with our code to put the wall as close as you like. The real and the imaginary part of the vacuum energy matrix V are defined by °°.
W~~1 = (fVX)’~”.
(7.19)
This expression must be compared with eq. (3.24)
W~’=
—~-
~
(a*i+h/2QjjaI+h/2i~0j+i/
(7.20)
2)~
and 2 @1~)1+hhI2= X,~1—
nq~-1”~i~_~I = HRX, 2 Eqs. (7.18)—(7.20) lead to VRI
=
X,
+
(aI)11h/
=
_______
—
L~x1
nq X,+ 1 +X1 2
=
H’X.
(7.21)
(7.22)
H*RIQHRI,
with —
~lJ
T
~xa+1/2
~xJ÷1/2
q2
n~ LWj+1/2
~.Wj+l/2
A~
Z4. Then—Oterms Similarly to (7.12), we expand (1/r’)v’ V’al’ into
~
1~112f
V’al’)
=~
d
1÷112,
(4w’. V’al’)
=~e1+h/2f;.+i12,
(7.24)
and obtain the matrix equations t~’Hd—Pe,
u=Rd—Se,
(7.25)
348
R. Gruber eta!.
/ ERA TO stability code
where t’ =
2
=
1
f
2sr
2ir
+
f~+i,2 27rL~0f dO
0
2ir
I
2irM
0
f f
HZ! = 2yr~0
2sr
21T
pit = ~
d01+112J —--v 0 ~dO,
f
d0’fj+112(4rr’G~5’),
f
d0’f,÷ii2(4rr’G~c’),
dO f~’+1I2 0 dOJ~+
2,r
1~2
2sr
2ir
f
~—~-~j 1 ~ dOf,+~p 0
dO’f,+112(4ir’G~1.),
0
2ir
I
v’. V’(4n~’G~5’),
2ir
2ir~O
=
~ r
2zr
f
dOf÷112
j
d0’f,+ii2(4~’G~~’).
Eliminating e in (7.25), we obtain 1R]1(r—PS~u). d [H— PS The vacuum energy (3.41) together with (7.19) becomes
(7.26) (7.27)
W~=
W~=
(7.28)
{X*(u+W+V)X]I,
with =
[~o
1~o1iJf~L~f] Sij, U
[d1dJ~o1~oJ/(~ d1~01)]o1J.
=
(7.29)
7.5. Resolution of the integrals The u integral in K~(~), eq. (7.1), can recursively be expressed in terms of elliptic integrals of the first the second (E) kind. It turns out that 4n ~+1 K0 2n+1 ~—1 2 K1=—K-i--—---—E, i~—l l—r~ —
K~+1
—
—
(K) and
2n—l K0_1, 2n+l
(7.30)
K0=K.
(7.31)
R. Gruber eta!. / ERA TO stability code
349
The n-derivative can be found by differentiating (7.24) and by using 1 ~K 1 ri + 3~ 0 1 ~E ~l E (3 + ii)Kj~ = K]. a~ 2Ø~ 1) i~ a~ 2(t~ 1) ~ —
——
—_[——
—
(7.32)
—
All the regular integrals in (7.14) are performed by means of a 4-point integration formula
f f
(7.33)
ab[f(—h~, o)+f(o, h~)÷f(h~,o)÷f(o,h~)],
f(x,y)dx dy
x—a y—b
4) which avoids any values
where = 2a~/i7~ and h~eliminating = 2b~./i7~. Eq. consists of an integration rule in 0(h on the h~ diagonal x = y, thus the (7.33) numerical singularities there. In the code, the singular part of 5afl~,eq. (7.11), is calculated analytically, i.e. 1 ~an~
=
~j+i
— —
01+1
~ a’~’~f
dOJ~÷
f
112 Ui
£
2
+
(7.34)
(0~+~O~)2]. —
dO’ 2 ln[b
Note that for the matrix B (7.14), p and p’ are on the plasma surface and at 0’ = 0, p’ = p a logarithmic singularity is obtained in (7.34). In the integral E, however, there is no singularity but for a wall near the plasma surface (when b is smaller than the mesh size) an apparent logarithmic singularity is built up numerically. For such a case it is advantageous to treat this integral analytically, to let us then solve 01+1
f Ui
o
0i+1
dO
f
do’ ln[b2
+
(0
—
01)2]
=
3(0
—
0
0~)2+ 4b(0
1+1
—
O~)arctan
Ui
2 in b2 + [—b2+ (O~+~ 0~)2]ln[b2 + (O~+~ O~)2]. (7.35) b Here b2 is considered to be constant in the interval O~~Z0 s~0~÷~ i.e. b2 = 0 for the matrix B (the value of the integral is then —(O~+~ O~)2[3 in(0 2 = (p~ Pc)2/PpPc atO = (O~+~ + 0~)/2for the matrixE. 1+1 singularities 0)2]) andinb (7.26). In the same way we treat the logarithmic —
—
+
,
—
—
—
—
7.6. Definition of the shell The position of the conducting wall is defined in subprogram SHELL. At the moment we measure the distance between the plasma surface and the shell by the quantity REXT-l, where REXT is given in the NAMELIST. The shape of the wall is defined by Pw =
Pp +
(REXT1)[p~(0)+ Pp(7T)l/2.
(7.36)
This means that we add a constant quantity in the radial direction to the plasma surface. It is evident that any other shape can be introduced by changing SHELL.,REXT > 1000 means that the conducting shell is at °° without the stabilizing line current along the toroidal axis.
8. ERATO 3 8.1. Organization ERATO 3 constructs the symmetric block structured matrices A and B which form the eigenvalue problem Ax = w2Bx. Only the upper symmetric part is stored as a one-dimensional table. One block contains all the contributions for a given s-interval. The blocks overlap due to the linear basis functions of X in the radial direction.
350
R. Gruberetal. /ERATO stability code
x,. 12 11
TT
48 47
x
46 45
x
33
44 43
x
1931
42 41
x
6~2B 40 l 15 27 39
x
2 1
l4~26 1325
x
Xl
V1V1
23 35
155
—
10 9 Xl Vi Vi X2 Y2 V2 X3 V3 ‘13 X4 V4 V4 X5 Xl Ti Bi I vi \\ vi X2 _________ B2 V2 V2 “ X3 __________ B3
________
21
8 l
I
6 5
________
1729
I
V3 X4
________ 156
1 3
0——-
“
14
B4
V4 X5
—
Fig. 6. Block structures of A and B for N~= 4.
38 37 X2
V2V2 X3
s V3V3
X4
Fig. 7. Global numbering in the
Y4V4 X5
(s, x) plane.
We obtain a matrix structure as shown in fig. 6.and the corresponding numbering of the unknowns is given in fig. 7. Two unknowns per mesh point are set for the real and the imaginary parts. In the case given we have chosen N~= 4 and N~= 5. One block has the size 8(N~+ 1) X 8(N~+ 1) which in our example corresponds to 48 X 48. Since only the upper half is stored 48 X 49/2 matrix elements define one block. The overlap goes over 2(N~+ 1) unknowns, giving 12 X 13/2 matrix elements. The length of the whole matrix which corresponds to the total number of unknowns is (3N 5 ÷1)(2N~+ 2). In the construction of the matrices we calculate all the possible contributions per mesh cell (16 X 16). Regularity conditions (only for the first block), periodicity conditions (for each block we eliminate all unknowns which lie in the domain ir <~< 2ir), and boundary conditions (for the last block only) are directly imposed on
XLj,l
~‘~‘‘/2.j.1 VL,ys ).1
4,..~ 3
2 1 ~LJ
12~16 11 15
‘1~~y2j
Fig.
8.
j’l
7
1014 9~l3 ‘~‘L.~/2,j
~i.l
0
X component
X
V and V components
0
physical components
6 ‘-‘5 XL,tj
Local numbering in one mesh cell.
R. Grubereta!. /ERATO stability code
351
these local 16 X 16 matrices. After that we add them up in a block. Two blocks are filled up consecutively. The overlapping part can then be directly added up. This is done in subroutine CONMAT. The blocks ofA are stored on disk file NDA and the blocks of B on NDB, see table 9. For the local 16 X 16 matrices containing all the contributions per mesh cell we introduce a new numbering system (see fig. 8). These numbers are related to the global ones (fig. 7) by means of the subroutine ~DEC). 8.2. Potential and kinetic energy The 16 X 16 matrices, i.e. the contribution to the potential energy matrix A per mesh cell, are built up in the subroutine (AHYBRD). In the integral (3.11)—(3 .13) we replace the vector components X, Y and V by values in a mesh cell, namely XRI
=
[X~’ + x~J÷, + X~,,~, + X~,~+11/4
axRh/ax = [X~J+
1+ X~’1,1+1 X~’ ax’~/as= [X,~ +X ,~+, X~’ —
—
—
—
2(X~+i Xi), —
X~1,1]/
x~/+i}I2(s
1+, se), —
yRI
=
VRI =
(Y~’~12,1 + Y~’~12,1÷1)/2, aYRh/ax = (Y~12j+1
—
Y~’112,,)/(x,+1
—
Xi),
(V)~12,1+ Vy~’~12,1+1)/2, av’~-’/ax=(V~’~12,1+1 V 1/2,J)/(xJ+1 —x,).
(8.1)
—
The positions of the nodal points in one mesh cell are explained in fig. 3 and the numbering of the unknowns in one mesh cell for the code in fig. 8. Note that eqs. (8.1) are identical to centred finite differences which will be applied to the integrals (3.12) of the variational form (3.11) for the mesh cell (s~~ S s1+1 ~ X~ ~
i’~ ~~-‘--+nX’~,
qax
qa~
ax’~ av” as ax I~= HXR /3XnqX’ ‘2 =
~,
jl
—
—
ax’ av’ ~,
as
ax
nq V’ ÷ i3x
~—
2 XR+a lnr2 VR ÷~
I~~a ln r as
ax
+
ax
k 2 3 4
5 7
~—,
I’~= HX’ + (3~nqXR+ nq VR + (3x
ax yR asnqY’
Table 9 IODSK3(k)
1
,
reads NEWRUN from NSAVE writes NAMELIST back on NSAVE reads vacuum contribution from NVAC writes matrix A on NDA writes matrix B on NDB reads equilibrium quantities
ax
~
as
ax
-~-
ax
+
~-,
as
(8.2)
R. Gruber et al. / ERA TO stability code
352 ,
a
a ln r2 a ln r2 ax v’+ ax
in r2
as
ax’ a v’ a “ ôs
ax
ax
i~’=x’~’. The quantity H = (q/T) a (T/q)/as
+
2j~O
2.The potential energy W,, (3.11) is then 5sr/l V~t’I
w~,= ~-ff~-dx [aI~2~J~2)
+
b(J~2~J~2~) + c(J~2+I’3~)+ d(I~2+I’)
—
e(I~2+Iç)].
(8.3)
In the same way one has to proceed for the kinetic energy, where no interactions between real and imaginary parts take place. The kinetic energy becomes
K = +ff~dx{f(XR2+X12) +g[(VR
+ yR
—
P~~)2 + (V1 + Y’ + flxx’)21 + h(YR2
+
~2)}
(8.4)
The variation is afterwards performed after all the nodal values. 8.3. Symmet7y conditions As mentioned before we only treat half the plane, and add half intervals in the symmetric half plane in order to be able to impose the symmetry conditions (3.18) which imply X’~(x)= —X’~(2ir x), yR(x) Y’~(2ir x), —
—
X’(x) = X’(2ir Y1(~)= —Y’(2ir
—
—
x)~
V’~(x)=V’~(2~r—x), V1(x)=—V’(2~i—x).
(8.5)
To fulfil (8.5) we introduce new variables XR( 2~ _X)rX’~(2iT—x)
+X’~(x)r0, x) V’~(x)= 0,
~‘(2n~ —x)
=x’(2ir —X)—X’(X) 0, X’(2z x) V’(2ir x) + V’(x) = 0,
VR(2~r x) VR(27r t(2ir—~)=Y’~(2ir—~)— Y’~(~)=0, ~“(2ir—~) Y’(2ir—~)+Y’(X)~r0, ?‘ which have to vanish. We first write eqs. (8.6) in matrix form
(8.6)
~Ux,
(8.7)
—
—
—
—
—
(x~X’~,X1, VR, V’, y’~,Y’)
and perform a matrix transformation on the eigenvalue problem Ax = w2 Bx which becomes (8.8) with
A
=
(U_i)T AlT1,
B = (U_l)TBIrS.
(8.9)
In the transformation (8.9) the block structures as well as the symmetries ofA and B are the same for A and ~&. In the code, we impose symmetry conditions cell after cell in subroutines (INTEGR) and
R. Gruber et al.
/ ERA TO stability code
353
& 4. Regularity and boundary conditions The regularity conditions (3.16) are easily imposed by X~=X’01=0,
(8.10)
0~j~N,<
throughout the subroutine (AWAY). In our code we have 2 different boundary conditions. The simple condition for a wall straight on the plasma surface (REXT < 1.0) becomes X1=X~1=0,
~
(8.11)
The other possibility is to prescribe a conducting shell at some finite distance (1.0
big value of REXT (REXT> 1000) simulates a shell at infinity. For these cases we add the vacuum energy matrix V, defined in (7.22) to the last block of the matrix A. This is done in subroutine (ADDVAC).
9. ERATO 4 9.1. Eigenvalue problem ERATO 4 solves an eigenvalue problem 2Bx (9.1) Ax’~c~., for symmetric, block structured matrices A and B. In addition, B is positive definite. The eigenvalue w2 of the problem can be positive (corresponding to an oscillatory, stable displacement), zero (the marginal situation) or Table 10 Parameters for VEKIT AL(3) AL(1) = ALO AL(2) AL(3) EPSCON EPSMAC NPIN(9) NPIN(1) = 0 NPIN(2) = N~ NPIN(3) = (3N NPIN(4) = NPIN(5) = NPIN(6) = NPIN(7) NPIN(8) = NPIN(9) =
5 + 1)(2N~+ 2) NITMAX NPRNT NDA NDB NDLT NDS
I 0 0
I I I
shifting
eigenvalue eigenvalue by normalization eigenvalue by Rayleigh quotient c.~,2 = (xT, accuracy of eigenvector machine accuracy
I
vector iteration required number of blocks
I I I I I I I
number of vector components maximum number of iteration steps printer output channel channel number for matrix A channel number for matrix B scratch unit number scratch unit and storage ofx
0 0 0 0
number of eigenvalues less than ALO number of iterations used number of converged components 0 all subdeterminants (A) regular = —1 subdeterminant (A) singular
NPOUT(4) NPOUT(1) = NPOUT(2) NPOUT(3) = NPOUT(4)
NEG NIT NCONV NSING
The eigenvectorx found at the end is stored on unit NDS and will be reread by ERATO 5.
Ax)/(xT, Bx)
354
R. Gruberetal. /ERATO stability code
negative (a purely damped or a purely growing, unstable mode). The eigenfunction x contains the displacement components arranged in the way described in fig. 7. The problem (9.1) is solved by an inverse vector iteration, namely by VEKIT, a subprogram of the block matrix library HYMNIABLOCK [10]. To have access to all the eigenvalues we perform a spectral shift. Instead of solving (9.1) we solve 2Bx, (9.2) Ax~ where A = A w~Band ~2 = ~2 ~4. By an adequate choice of w~(the NAMELIST quantity ALO) we are able to pick out each eigenvalue of the spectrum that we like. To make the iteration efficient we decompose A into —
—
A=LDLT,
(9.3)
where L is a left hand side trigonal matrix with ones on the diagonal and D is a diagonal matrix. The number of negative terms in D is identical to the number of negative eigenvalues. Therefore we always know which eigenvalue of the spectrum has been calculated.
9.2. How to call VEKIT The eigenvalue solver VEKIT is called in EIGVAL by CALL VEKIT(AL,EPSCON,EPSMAC ,NPIN ,NPOUT).
(9.4)
The meaning of the formal parameters are described in table 10 (I 10.
=
input parameter, 0
=
output parameter).
ERATO 5
10.1. Diagnostics This main program performs all the output with the exception of the eigenvalue which is printed out in VEKIT, a subprogram of the library HYMNIABLOCK [10]. We first read in the eigenfunction from disk unit NDS, separate it into the six components XR, x’, yR y’ VR and 0 which are transformed back into ~ and see eq. (3.9). These components are normalized with its highest absolute values (NORM FACTOR * 999) and the three most significant digits are printed out as a table ins and x. The s-coordinate goes from up (s = 0) to down (s = 1) and the x coordinate from the left (x = 0) to the right (x = ir). The actual values of the vector components are obtained by multiplying the tables by the corresponding NORM FACTO1~which is printed out at the top of the tables. As a next step we plug the eigenmode back into the variational form, eqs. (3.l1)_(3.15), and calculate the contribution of each mesh cell. Again we print out a table containing these cell contributions to the potential energy and to the kinetic energy. The last column at the right hand limit contains the cross sum of all x contributions for a given s interval. With such a diagnostic we are able to see where the negative driving terms of an unstable displacement come from. The test case in section 12 shows such a plot for an internal kink mode (all the negative contributions to the potential energy are internal). As an additional control, we add up the potential and the vacuum energy and form together with the total kinetic energy the Rayleigh quotient which should correspond to the eigenvalue obtained by HYMNIABLOCK. The diagnostic part concludes with a “MAP OF POTENTIAL ENERGY CONTRIBUTIONS” which enables us to find out which of the 5 integrals I~to I~contributes to the potential energy. ~,
~,
~,
~
~,
10.2. Plotting The plotter or graphics output is controlled by the quantity NFIG (in NAMELIST). NFIG = 0 no plotter output, NFIG ~ 16 displacement graphs as shown in the test run output are drawn at NFIG different angles at which the torus is cut perpendicularly to the magnetic axis (in the s—x plane),
R. Gruber et al.
/ ERA TO stability code
355
ARROW enables the size of the arrows to be changed, ANGLE (NFIG) denotes the angles in degrees at which the NFIG torus cuts have to be drawn. ANGLE is given in NAMELIST. The quantity ITEST is CDC dependent. It must be set to 1. ITEST = 0 only checks if the plotter output is in the limits without drawing anything. The plotter or graphics output passes through the following subroutines ABSDIR calculates the directions and the amplitudes of the displacement in each mesh cell, PLOTER organizes the plot, FLECHE draws an arrow, SURF draws the plasma surface, TITRAX produces titles and axes. The plotter or graphics output contains three machine dependent subroutines which must be adapted for each computer: PLOTS initializes the drawings, PLOT moves the pen up or down, SYMBOL draws a string of characters. Details about these three subroutines are given in the listing of ERATO 5. To avoid using another machine dependent subroutine (NUMBER) digits are converted into characters by ENCODE which is not standard FORTRAN but available on many computers. The character obtained is then drawn by SYMBOL. The CDC subroutine DATE is also used.
11. Test cases 11.1. Preset values (internal mode) The code can be run without giving any data between the $NEWRUN and the $END card. This default case calculates an internal mode of the Solovev equilibrium, eq. (2.8), with the parameters given in table 5 of section 6. As a result one obtains an eigenvalue w2 = —1.5883. Note that the $NEWRUN card must be preceded by 4 label cards. Setting NTCASE = 1 in the namelist gives the result w2 = —1.5889.
11.2. Test run output (kink mode) As a second test case we treat a free boundary (REXT = °°)kink mode with a growth rate _~2 = 1.2373. This case is shown as test run output at the end of the paper. Only the first two and the last two pages are shown. As input card stream the card sequence 1st card: this is the ERATO test case 2, 2nd card: an external kink mode (n = 2, E 2 is calculated), 3rd card: note that these 4 cards must always 4th card: precede the NEWRUN input, 5th card: $NEWRUN, 6th card: NPSI = 14, NCHI = 15, WNTORE = —2., 7th card: ELLIPT 4., NTCASE 2, ALO —1.2, 8th card: QIAXE 1.6667,REXT= 1001.0,NMESH= 1, 9th card: NFIG = 4, ANGLE(1) 0., 30., 60., 90., 10th card: $END must be given to obtain the test run output. For this case we also show the plotter output of the resulting typical kink-like eigenvector (fig. 9). The (~,~) displacements are shown at each (s, x) mesh cell centre at (NFIG = 4) differentq angles (ANGLE(l)= 00, 300, 60°,90°).Changing REXT = 2.0 gives an eigenvalue w2 = 1.1015. I~2
356
R. Gruber eta!.
~ /1 ~ IT~V.’,f~~
code
/~\‘-/~\
ii
/
--
/ ERA TO stability
‘
~ IV
/‘~,,/:
~
/~st~
~
~
~
~
:
~‘
~ço ~
:::::: --
~
~:-
~ I
~
LPUSPNNE
-
/L
~
~
V
if
/
- -
~~
“,
I
-
L6US6NNE
.
S
‘
,
—
.•‘~
:~O~j , c~—~:~_
5
\‘
— —
—.~
~
\~ ~H~~•’’’-\L’ ~
1
\V~~r
/
-
/
/
/
f
~EP6TO
=
/
/
~j~=-
60
26/03/86
/
-
/
\V
P1-/I
—
-
_\
\
~
rpr
26/03/80
—\
——S
// / / 1y~~//~ ~ / V \l~ .-‘~1~’/ ‘I ,‘~L~1~1 ‘I \~ S~,y
~l4
30
-N-\
~
\
/1~
I
‘ERPTO
~
~~‘-~‘‘
I
LRUSRN~IE
—
~ /~ ~~ I’H~~\
V I~~~V_: ~
CPPP
-~
IV
_—_/7~
~
\V
EPF
CPPP
~\
~
-‘ -
-
26/03/80
/
-
~‘-.-,,
1
V
/ -
:~~
PHI
‘S
L’-
H
ERPTU~
~
\
=
‘
7.I~
\
‘~ ft
I V~’~ ~
£PF
,/
~“ -
‘~
_
p~
CRPP
~
-
:~ ‘
\~
~
CRPP
£PF
—
LRUSRNNE
PHt
‘cR610
Fig. 9. Plotter output showing the displacement vector field at different ~ angles.
=
90
26/03/86
R. Gruber et al.
/ ERA TO stability code
357
All these cases are related to the calculations performed in ref. [27]. The computing time of such a case is around 2 mm on our Cyber 170-720.
12. Instructions to the user
12.1. Convergence properties A regular finite element approach [20] always shows quadratic convergence behaviour when the number of intervals goes to infinity. As presented in a previous publication [8] we have to take a very large number of mesh cells to obtain the required quadratic convergence and, sometimes, the actual computers are not big enough to solve the stability problem by such a method. This inconvenience comes from the fact that the eigenvalue is always approached from the stable side of the spectrum, i.e. from above. A coupling between the unstable, not yet enough converged mode, with stable, almost marginal modes of the continuous spectrum, can strongly affect the convergence behaviour. In the finite hybrid element method, however, we reach a lower energy level, and therefore prevent the coupling with the continuum. But another problem arises which we call the spectral shift in q [15]. However, the elimination of this spectral shift, eq. (4.17), is only valid in the case of a Solovev equilibrium. We have seen that this shift depends on the shear at the plasma surface. A more general expression has still to be found. In fixed boundary cases we even sometimes find higher order convergence laws. Such a higher order convergence is only found in the absence of the vacuum region since the vacuum contribution in its actual form always shows a quadratic convergence law. For Tokamaks it is seen that the maximum F2’s are always higher in the discrete case than for the exact value. The marginal situation of the continuous spectrum can also become slightly unstable and a convergence study must push this apparent unstable mode back to the stable region. A careful study of the localization of such a mode by means of the potential energy and the kinetic energy tables can be a helpful tool to recognise a continuous mode in the unstable region. As a consequence, the user must be very careful in the interpretation of his results. 12.2. How to switch from equilibrium to ERA TO As mentioned before, ERATO needs to know t~(r,z), p(i,li), T(’I’) and TT’(i~),such that the equilibrium equation = —r2p’ TT’ (12.1) —
is fulfilled. Whenever the equilibrium code is not written in this “natural” system [18], one has to perform a transformation of p’ and TT’ in the subprogram TANDP. The normalization with respect to the magnetic axis R and the normalization of T with respect to the magnetic flux at the axis T 0 is done in subprogram NOREPT -and the user has nothing to care about. As a default value, we use unit NDORY = 11 (catalogued on local file name DORY) to pass the necessary quantities from the equilibrium code to ERATO. The accuracy of a numerical equilibrium using a finite difference method is somehow related to the accuracy of a finite element method as used in ERVATO. As a rule of thumb this relation is essentially given by Nr
6N5,
N2
3M5.
(12.2)
These relations come from the fact that in a finite difference approach the derivatives are found in between two intervals, whereas in the finite element method only one interval is used. In addition the number of r-intervals must be doubled in respect to the number of z-intervals, since only the upper part of the torus is treated. Another
358
R. Gruber et a!.
/ ERA TO stabi!ity
code
factor of 3/2 is introduced to take into account that on the average about 1/3 of all the r—z mesh cells are outside the plasma region. Note that in the high i3 case where the magnetic axis is pushed outwards the factors in (12.2) must be increased. 12.3. Choice of NAMEL1ST parameters All the NAMELIST parameters can be modified in all the 5 main programs. Generally one defines them in ERATO 1 and leaves them unchanged throughout the whole tr~atment.However, it can sometimes be of some interest to build the matrices A and B only once and go through the spectrum with different values of ALO c4. The matrices A and B are then stored on disk files NDA = 7 (with local file name TA) and NDB = 9 (TB) and the NAMELIST quantities are catalogued on unit NSAVE = 8 with the local file name OUT. In the NAMELIST input, there are some special points to be considered: (I) Convergence studies are simpler when the same number of intervals in s and x is chosen, i.e. when N~=
N5
+ 1. (II) For an unstable configuration ALO <0. (III) A fixed boundary case yields REXT < 1; REXT> 1000 means a wall at infinity and in between there is a wall at a distance (REXT 1) X (plasma diameter at z = 0)12. (IV) The quantity EPSMAC has to be set to the machine accuracy. For a CDC 6500, EPSMAC = 10_12. (V) The mass density p is preset to i in the whole plasma region. When a mass density profile p(s) is required one has to foresee this in the subprograms EQANAL and NUMEQU when filling up EQ(7,JCHI). The user —
EPF LAUSANNE MOS/BE 1.3 1EV 48~—0621f07180 12.12.28.72561gM FROM 12.12.28.IP 00071232 WORDS — FILE INPUT DO QL 12.12.28.T2551,CMILi0000,Ti000. GRUBER,I~RPP 12.12.28. 1 12.I2.30..REQUEST,SOUROE,~PF. 12.12.30.UPOATE.F.N=OLDPL,C=0,L=0,S. P i2.12.30.REPARES OLDPL 2 12.12. 31..UPDATE OREATION RUN 12.i2.3i.CREATIN~ NEW PROGRAM LIBRARY 12.12.Li9. UPDATE COMPLETE. 12.12.49.CATALOG, SOURCE,ERATO, IO=LISRS. PP=999. 12.12.’.3.INITIAL CATALOG 12.12.51.CT 10= LIARS PFN=ERATO 12.12.51.CT CY= 001 00020 PBS 12.12.5I.UPOATE,Q,L=Q. R 12.12.51.OUTINES FOP ERATOI 3 l2.~3.10. UPDATE COMPLETE. I2.~I3.I0.FTN,I=COMPILE,L=0. 12.13.10. Li 12.i’..kS. 32.587 OP SECONDS COMPILATION TIME 12.14.49.MAP(PART) 12.1~.Li5. 5 12.14.le.LGO. ~ ERATO I ~ 12.1Li.1.5. S 12.15.22. STOP 12.15.22.0554000 CHAMP MEMOIRE MAXIMUM 12.15.22. 27.007 OP SECONDS EXECUTION TIME 12.j5.22.RETURN,LGO. 12.15.22. 7 12.15.22.UPOATE.9,L0, R 12.i5.22.OUTINES FOR ERATO2 8 12.15.28. UPDATE COMPLETE. 12.15.28.FTN,ICOMPILE,L0. 12.15.28. 9 12.16.08. 12.095 OP SECONDS COMPILATION TIME 12.I6.08.LGO. ~~4” ERATD 2 12.16.08. 10 12.16.13. STOP 12.16.13.03560/lB CHAMP MEMOIRE MAXIMUM 12.16.13. .145 OP SECONDS EXECUTION TIME 12.16.1 3.REWIND,LGO. 12.16.13. ii i2.i6.13.UPOATE.9,L0. P 12.15.13.OUTINES FOR ERATOS 12 12.16.19. UPDATE COMPLETE. 12.16.19.FTN,I=COMPILE,L=0. 12.16.19. 13 12.16.54. 11.887 CP SECONDS COMPILATION TIME 12.16.54.LGO. ERATO 3 12.16.54. jL 12.17.22. STOP Fig.
12.17.22.OT6k008 CHAMP MEMOIRE MAXIMUM 12.17.22. 10.880 CP SECONDS EXECUTION TIME 12.17.22.REWINO,LGO. 12.17.22. 15 12.t7.22.UPOATE,Q,L0. P 12.i7.22.OIJTINES FOR ERATO4 16 12.17.29. uPDATE COMPLETE. j2.17.29.FTN,I=COMPILE,1=0. 12.17.29. 17 12.18.06. 6.362 OP SECONDS COMPILATION TIME 12.18.06.FILE.TAPEIO.RTS. 12. 13.06. 18 12.18.06.LDSET.FILESTAPEIO. 12.18.06. 19 12.18.06.LGO. ~ ERATO L 12.18.06. 20 12.22.36. STOP 12.22.36. 86.5L+L. CP SECONDS EXECUTIOP~TIME 12.22.36.i0Lii000 CHAMP MEMOIRE MAXIMUM i2.22.36.REWIND.LGO. 12.22.36. 21 12.22.36.UPDATE,9 1=0. P 12.22.36.OUTINES ~0R ERATO5 22 12.22.Li7. UPDATE COMPLETE. 12.22.Li7.FTN. I=COMPILE,L=0. 12.22.47. 23 12.23.19. 14.252 OP SECONDS COMPILATION TIME tZ.23.l9.LOSET,FILESTAPEIO. 12.23.19. 24 12.23.19.IGO. (RAT/I 5 12.23.19. 25 12.23.24. NOil~FATALLOADER ERRORS — SEE MAP 12.23.32. STOP 12.23.32.0644008 CHAMP MEMOIRE MAXIMUM 12.23.32. 2.604 OP SECONDS EXECUTION TIME 12.23.32.EXIT. 12.23.32. 26 12.23.32.QP 00010816 WORDS — FILE OU1PUT, DC Lie 12.23.32.MS 86016 WORDS ( 526848 MAX USED) 12.23.32.CPA i84.q58 SEC. 6.165 AD.). 12.23.32.CPB 59.470 SEC. 1.982 ADJ. 12.23.32.10 160.905 SEC. 10.726 80.). 12.23.32.C$ 1075i..317 045. 388.972 AOJ. 12.23.32.SS~-~COIJT OU CALCUL 285.90 FR ~‘ 12.23.32.00 PERIPHERTE EN SUS (MESS CC. TRACEUR) 12.23.32.PP 308.976 SEC. DATE 08/09~80 12.23.32.CC 25.70 FR I/O 0 PAUSES 12.23.32.EJ END OF JOB. ~‘ rn~8~4*~#
10. Control stream of ERATO.
T2SSI9M 125519M
I//f If/f
END OF LIST I//I END OF LIST If/I
R. Gruber eta!.
(VI)
(VII) (VIII)
(IX) (X)
/ ERA TO stability code
359
is warned against too realistic mass density profiles which fall to zero at the plasma boundary. These can cause unrealistic values of w2 which tend to infinity with 1/p (see ref. [19]). NLY = T calculates the whole problem, NLY = F eliminates the toroidal component Y by imposing V• = 0 in the potential energy and Y = 0 in in the kinetic energy, NLTORE = T for toroidal geometry (tokamaks, spheromaks, reversed field pinches), NLTORE = F for helical geolnetry. This case is not yet included in the published version of the code. NMESH = 0 for an equidistant s-mesh, NMESH = I for a non-equidistant s-mesh. Accumulation around integer nq surfaces is made. The user can easily change the accumulation algorithm in MESH. WIDTH 00 for equidistant x-mesh, WIDTH = finite accumulates the x-mesh at the outer region. The user can again easily change this. 1—0 operations in IODSK3 of ERATO 3 and in IODSK4 of ERATO 4 should be accelerated by introducing direct access of the matrix block records and by directly buffering into main memory.
10.01
ERATO1
1.1.02
CLEAR
(cocJ LOCE
11.03
PRESET 1.1.01
LABRUN 1.1.04
DATA 1.1.05
AUXVAL 1.2,01
ORGNIZ 2.0.01
ERATO2 3.0.01
ERATO3 4.001
ERATO4 5.0.01
ERATO5
2.2.01
VACUUM 3.2 01
AANDB 4 2.01
EIGVAL 5.2.01
SORTIE
Fig. 11. Overall structure of ERATO. The main programs ERATO2—ERATO5 are dummies, and the subroutines VACUUM— SORTIE can all be included in ERATO1 if sufficient main store is available.
360
R. Gruber eta!.
/ ERA TO stability code
Table 11 Adjustable COMMON dimensions No.
Block name
Variables
Symbolic dimensions
2.6
COMVEC
XT, UT, VT
(6N
X,U,VA,VB,C A1-A6 BETRZ, CHIRZ RRO, RRS, RZS, TETA
NRZ
131
lNr,Nz] Nr
(65,33)
9.1
COMMAP
CPSEQ
EQPS,EQR,EQT,EQZ ILOC 9.2
COMVID
8(Nx+l) 5 + 2)(Nx
+
1)
Default a) 1376 128
65
15
Nx
G A—F WVAC H
[2Ivx
[4Nx_2,4Nx_2] (2Nx 2)(2Nx 2)
— 2, 2Nx — (2~’~T~ — 2, 2Nx —
—
2,61 2]
—
9.3
COMMTR
A, B
4*(Nx + l)*(8Nx + 9) + 1
9.5
COMPLO
CNRRZR X XS,YS KK~KW~RK~RW}
ENs+1,Nx+l]
DW
[Ns,Nx+l,61
SUM
[N
(6Ns+2)*(Nx+l) 2*NRZ ENs+1,Nx+1J
5+1,6}
(28, 28, 6) (28, 28) (58,58)
784 8257 (15,16)
1376 262 (15,16) (14,16,6) (15,6)
a)Dej~auItforN~.:65,N~=33N~ l4,N~~’ 15,NRZ= 131.
12.4. Control stream and data set As an example (fig. 10) we show the control stream of the Solovev test case. All the 5 ERATO main programs are corrected by a batch editor UPDATE, compiled and executed. The control cards preceding ERATO 4 and ERATO 5, i.e. the FILE and the LDSET cards are set to eliminate problems with the BACKSPACE command in IODSK4. After the control stream, the UPDATE modifications for all the 5 main programs follow between the first and the second correction decks the NAMELIST preceded by the 4 label cards is inserted. They are read in ERATO 1. All these card decks are separated by end-of-record cards.
12.5. Dimensions The dimensions which have to be arranged are given in table 11. In the original version the default values are = 14,N~= 15,J\1r65 andN5 =33. These arrays have also to be modified when bigger cases should be treated. For the surface fit we have foreseen NRZ ~60. N5
R. Gruber et aL
/ ERA TO stability code
361
Table 12 ERATO 1, subprogram names and identification Name
No.
Number of
Contents
arguments
LABRUN CLEAR PRESET DATA A(JXVAL ORGNIZ FORMEQ MESH KERTP KERTES SCALES FITAXE NOREPT TANDP SELECT ISTPSI EQANAL ALPHA
1.1.01 1.1.02 1.1.03 1.1.04 1.1.05 1.2.01 1.2.02 1.2.03 1.2.04 1.2.05 1.2.06 1.2.07 1.2.08 1.2.09 1.2.10 1.2.11 1.2.12 1.2.13
PSRHO
1.2.14
5
labels the run clears variables and arrays set default values defines data specific to run set auxiliary values organizes the calculation organizes equilibrium mapping set up mesh prepares for NTCASE = 1 (Solovev analytic) prepares for NTCASE = 2 (Solovev numerical) scales the numerical equilibrium fits equilibrium around axis normalization to ERATO standard set T(~~), p(9) and p(p) selects analytic or numerical equilibrium radial counter of (v,, z = 0) analytic equilibrium quantities calculates q and non-orthogonality givesrandz(s, 8) ~p(r,z) by fit r, z derivatives of 1~/(r,z) interfaces numerical equilibrium — ERATO advances on s = constant line interpolates on s = constant
1
set up numerical EQ
1
PSIRZ
1.2.15
1 2 1 2 5 2
DERIVE
1.2.16
2
INTFAC
1.2.17
1
ADVANC INTPOL NUMEQU
1.2.18 1.2.19 1.2.20
5
SURFIT TRIBAC MERCIR INTGRA QUADRA
1.2.2 1 1.2.23 1.2.24 1.2.25 1.2.26
BALOON ALDLT
1.2.27 1.2.28
5
BETAST
1.2.29
1
DOMAIN SMLEQU
1.2.30 1.2.31
YYP BNDRY
1.2.32 1.2.33
7
OUTPUT IODSK1
1.3.01 1.3.02
1 1
6 6 3
s, x
fits the surface solves Ax = b Mercier’s criterion calculates integrals for Mercier criterion polynomial fit for MERCIR
ballooning criterion decomposes A = LDLT beta values domain for analytic fit
analytic equilibrium in domain parabolic fit dummy routine prints output disk operations for ERATO 1
12.6. Space and time requirements The memory space of ERATO can be estimated by the formula
W~l6OOO+(N~+l)(32iVx+ 18N8+42).
(12.3)
This yields 49000 60-bit words on a CDC 6500 for a 24 X 25 case. The computing time is given by T~64N~Ns(tiN~ + t2N1~+
t3),
(12.4)
R. Gruber eta!. / ERA TO stability code
362
11213
~LJ~ J1224
2O3
Jl
I
1225
1~.1~C!R
MESH
INTGRA
Jl.2.04
I
Ji
1210
I
SELECT
11302
]
7
J1.2 05 I KERTES
302
1201
ORGNIZ
1252 FORMEQ
1219
IODSKI
INTPOL
]
1 3.02
BND~J
IODSK1
1D~IVE
1
ISTPSI
J1218
7
I
Ef7
1
PSRHO 11216
I
11215
fl~7PSIRZ
7 DERIVE
1 2 29
I 7
~‘l
ADVAI’C
[1214
BETAST
11302
I 7
IODSIO
II 223 7~TR1BAC I 11.219
7
INTPOL
1211
DOMAIN
ISTPSI
1207
7
IODSK1
H
VYP
302
1,233
J1230
__________~
232
NUMEQU
ALPHA 1126
QUADRA
7
220
7
~
[i~211
[j~~J~
KERTP
7
1226
~1
I
1u~~
T 1
M3REPT
l12~
12.15
7 SMLEQU
_________
J1206 ~
I
SCALES
_Ji
2.06
I_________ 112
21
~‘[
LSURFIT
J7
]12.17 ~ INTFAC
]13 01
27
L
112 28 ALDLT
1223 TRIBAC 1219
BETAST
_________
Fig.
where ~
ISTPSI
ADVANC
11229
7
211
1 2.18
7 OUTPUT
7 SALOON
Ji
12. Subroutine map of ERATO1.
is the number of iterations in the eigenvalue solver. The time constants are t
r2
1
23 .ts when machine
coding is used for the matrix decomposition, 50 jts and t3 150 is on a CDC 6500. For a iest case with 24 X 25 intervals, the decomposition takes 9 mm, the 4 iteration steps take 3~min. All the rest (i.e. mapping,
vacuum, matrix construction, shift and diagnostics) take 2~-min. Altogether, such a case takes 15 min on a CDC 6500 which corresponds to about 1 mm cpu time on a CDC 7600.
R. Gruber eta!. / ERA TO stability code 2.3.01
363
MOOS
IODSK2 2.2.02
CALD j22O4
TETMSH 2.2.03
2.2.06
ABCDEF
[__ji~
ROTETA {
CONCEL
SHELL
] r~i
2.2.12
EKN
COC
[2.2.0) VACUUM
LOCF 2.2 07
2.2.08
QCON 2.2.09
2.2.14
WCON 2.210
IMOC 2.2.11
SYMETR Fig.
ZERQJ
MULT
13. Subroutine map of ERATO2.
Table 13 ERATO 2, subprogram names and identification
Name
No.
Number of
Contents
arguments VACUUM TETMSH ABCDEF ROTETA SHELL CONCEL QCON ZERO
2.2.01 2.2.02 2.2.03 2.2.04 2.2.05 2.2.06 2.2.07 2.2.08
WCON
2.2.09
SYMETR
2.2.10
MULT
2.2.11
4
product of two matrices
EKN EK IMGC LODSK2
2.2.12 2.2.13 2.2.14 2.3.0 1
4 3 4 1
elliptic integrals completes elliptic integrals inverts a matrix
organizes vacuum calculation
3 5 3 1
8-mesh in vacuum set up matrices it, B, C, D, E, F prepares p. 0 fdl~integration defines position of the shell contribution per cell
calculates Q-matrix reduces for n = 0 total vacuum energy matrix imposes symmetry conditions
disk operations for ERATO 2
[
1
EK
R. Gruber eta!. / ERA TO
364
stability code 3.2.10
DEC 3.2.04
AHYBRD 3.3.01
3.3.01
1005K3
3.2.01
3.2.02
AANDB
100SK3
3.2.03
CCt’JMAT
INTEGR
3.2.05
BHYBRD 3.2.07
TRANSF 3.2.09
AWAY 3.2.06
STORE 3.2.08
A~VAC Fig.
33.01
IODSK3
14. Subroutine map of ERATO3.
Table 14 ERATO 3, subprogram names and identification
Name
No.
AANDB CONMAT INTEGR AHYBRD BHYBRD STORE
3.2.01 3.2.02 3.2.03 3.2.04 3.2.05 3.2.06
TRANSF
3.2.07
ADDVAC DEC IODSK3
3.2.09 3.2.10 3.3.0 1
Number of arguments
3 2 2 2 5 3 3 1
Contents
organizes matrices A and B creates blocks of A andB cell contribution and conditions 16 X 16 function forA 16 X 16 function for B adds cell contribution transformation ofA and B removes a column and a row address correspondence disk operations for ERATO 3
Table 15 ERATO 4, subprogram naitses and identification Name
No.
EIGVAL SET3
4.2.01 4.2.02
Number of arguments
Contents
organizes eigenvalue calculation set quantities for VEKIT
3.2.09
AWAY
R. Gruber eta!. / ERA TO stability code
365
MOO?
FXFU MDO3
M012
IODSK4
BSI-UFT MDO4
MOOS
CONALR
CALD
MD12
MD12
IODSK4 ~
4.2.02
MDO6
5ET3 j4.2.O1
IODSK4 MOb
CBXLU
GETRG
MOO)
EIGVAL
M011
VDUT
PUTRO M007
CDLHXV
M012 .
IODSK4 MO1O
GETRG M008
MOll
CONORM M009
PUTRG M012
EIGEN
IODSK4 MOlO
GETRG Fig. 15. Subroutine map of ERATO4. Table 16 ERATO 5, subprogram names and identification Name
No.
SORTIE DEBUT BACKTR ABSDIR VACPOT VALUES DIAGNO PLOTER
5.2.01 5.2.02 5.2.03 5.2.04 5.2.06 5.3.01 5.3.02 5.3.03
FLECHE SURF
5.3.04 5.3.05
TITRAX
5.3.06
Number of arguments
Contents
organizes printing and plotting
prepares vector components back to real displacements prepares arrows
1
3
calculates vacuum contribution prints value table diagnoses energies organizes plot plots an arrow plots the plasma surface plots title and axes
366
R.
Gruber eta!. / ERA TO stability code 5.2.02
3.2.10
DEBUT 5.2.03
BACKTR 5.3.01
VALUES 5.2.01
5.2.04
SORTIE
ABSDIR 5.3.03
PLOTER
DEC 7
I
(GRAPHICS] .
PLOTS [GRAPHICS]
PLOT 5.3 04
FLECHE 5.3.05
SURF 5.3.06
TITRAX 5.3.02
DIAGNO
5 2.06
VACPOT
~RAPHICS]
PLOT (GRAPHICS]
PLOT (CDC]
DATE [GRAPHICS]
SYMBOL
Fig. 16. Subroutine map of ERATO5.
Acknowledgements We would like to deeply acknowledge the continuous support by the applied mathematics group at EPF Lausanne, especially Prof. J. Descloux, Drs. H. Evéquoz, Y. Jaccard and J. Rappaz. Fruitful discussions with Drs. T. Takeda and T. Tsunematsu brought one author (R.G.) to the correct formulation of the finite hybrid elements. A strong competition with the Princeton MHD group forced us to do our best. Many helpful propositions and critics were given by Drs. K. Appert,~L.A.Charlton, R.A. Dory and J.L. Soulé. All this work would not have been possible without the excellent service given by the Computing Center at EPF Lausanne. This work was supported by the Ecole Polytechnique Fédérale de Lausanne, The Swiss National Science Foundation and by Euratom.
References [1] A. Sykes and J. Wesson, Nucl. Fusion 14 (1974) 645. [2] G. Bateman, W. Schneider and W. Grossman, Nuci. Fusion 14 (1974) 669. [3] M.S. Chance, J.M. Greene, R.C. Grimm, J.L. Johnson, J. Manickam, W. Kerner, D. Berger, L.C. Bernard, R. Gruber and F. Troyon, J. Comput. Phys. 28 (1978) 1. [4] I.B. Bernstein, E.A. Frieman, M.O. Kruskal and R.M. Kuisrud, Proc. Roy. Soc. (London) Ser. A244 (1968). [5] K. Appert, D. Berger, R. Gruber and J. Rappaz, J. Comput. Phys. 18 (1975) 284. J. Rappaz, Numer. Math. 28 (1975) 15. [6] R.C. Grimm, J.M. Greene and J.L. Johnson, in: Method of computational physics, vol. 16 (Academic Press, New York, 1975) ch. 4. [71 R. Gruber, 1. Comput. Phys. 26 (1978) 379. [8] 0. Berger, R. Gruher and F. Troyon, Comput. Phys. Commun. 11(1976) 313. [9] L.S. Solovev, Zh. Eksp. Teor. Fix. 53 (1967) 626; JETP 26 (1968) 400.
R. Gruber et al. / ERA TO stability code
367
[10] R. Gruber, Comput. Phys. Commun. 20 (1980) 421. [111 J.D. Callen and R.A. Dory, Phys. Fluids 15 (1972) 1523. [12] B.B. Kadomtsev, Reviews of Plasma Physics, Vol. 2 (Consultants Bureau, New York, 1966) p. 155. 1131 E. Rebhan and A. Salat, Nuci. Fusion 16 (1976) 805. [14] K. Appert, D. Berger, R. Gruber, F. Troyon and K.V. Roberts, Comput. Phys. Commun. 10 (1975) 11. [15] R. Gruber and F. Troyon, Proc. IRIA, Paris, eds. R. Glowinski and I.L. Lions (Springer, Heidelberg, 1979); Lecture Notes Phys. 91(1979) 288. [16] M.H. Hughes, K.V. Roberts and G.G. Lister, Comput. Phys. Commun. 10 (1975) 167. [17] D. Dobrott, D.B. Nelson, J.M. Greene, A.H. Glasser, M.S. Chance and E.A. Frieman, Phys. Rev. Lett. 39 (1977) 943. J.W. Connor, R.J. Hastie and J.B. Taylor, Phys. Rev. Lett. 40 (1978) 396. [18] E.S. Weibel, Am. J. Phys. 36 (12) (1968). [19] K. Appert, D. Berger and R. Gruber, Phys. Lett. 46A (1974) 339. [20] O.C. Zienkiewicz, The finite element method, 3rd ed. (McGraw Hill, London, 1977). [21] S. Rousset, R. Gruber and F. Troyon, ZAMP 30(1979) 862. [221 F. Troyon, L.C. Bernard and R. Gruber,Comput. Phys. Commun. 19 (1980) 161. [23] R.G. Bateman and Y.-K.M. Peng, Phys. Rev. Lett. 38 (1977) 829. [24] K. Lackner, private communication. [25] T. Takeda and T. Tsunematsu, to be published. [26] C. Mercier and H. Luc, The MHD approach to the problem ofplasma confinement in closed magnetic configurations, Lectures in Plasma Physics (1974), EUR 5127e. [27] D. Berger, L.C. Bernard, R. Gruber and F. Troyon, ZAMP 31(1980)113.
368
R. Gruber eta!.
/ ERA TO stability
code
TEST RUN OUTPUT THIS IS rr1+ AN
k 4 1
3
I9~,+~2I
IXTISNAL KIN~( Iluti
NOTE
TruSt
THAT
PRECtuE Tic —.17654*06
H
TIlT CASo IS LALCILATIU
4
LAI2S
lUST AL,AYS
A
S
1 S I
I
iNPuT
_.6U9A1E~Li —.104791*0J
—..,?01r,t—ruI .I30354*i0
—.3545c+—0Z —.74s33u.—.~l. .112L’44
—.LluO{)E*40
—.c132?t*u4
—.17504,..Uu
—.184881*LrU —.I54941*uO —.69tj~4E—rJI
—.Id4A4L+.r~ —. 9Ji*uu —.bi*,i—uI
—.10113.040 —. I ISocIc —.3+4101—41
4
26
.3,7b41.~r.r2
31
—.10940+—ui —.1ru1,6L—UI —.213531—0+—. 07905+—ui .L bb+’u., .141041+lu —.I7043~~u* —.isu~,.cuu —. *.riL2o...JJ —.+7*91.uu —.i391l 1ou—.~. 27u,I+oj —. 23531+—ni —.013751—nil
—.254t~.L—r.r~. —.944601—ui —.lJrr*jri —.i021c*uu —.I7042+*u0 —.II544c~uu .o1415c—U2
—.326i71—uL —.1or3~3E*0u —.IS~r.u01*uJU —.1d335**UU —.473~3L*44 —.Iu9,+5+*Uu
—.397731—UI —.lU7Ibl*OU —.159221*00 —.i84301+00
—.15950**r30 —.S9LSbE—01
—.468811—lA —.11313c*UU —.163161,00 —.185U51*urJ —.165221*06 —.884611—Ui
~.53930lr01 —.II9U5E~uU —.1.66831*04 —.L83i.8E~U0 —.16037E*U*r —.770191—01
is
40
FIr
I
2
nn.
4. i84Eu1
9.8241—11
in.
no.562+—uI
0.5101—ui. 7.6431—UI 3.7754—UI
b.b24E—OA 7.7561—01. 8.8891—61
9.793+—Li
3.9081—UI
L.uOZL*UU
i.i33c+nnnn 1. 700+ *00
+.1041*U0 £ • 2171 *~JU
1.4151*00 1.2291 *40
5.7371—UI 7. 869s—UI 9.0021n—U1 1.0131*00 1.1271*60 1. 2401 *00
4.1791—LU
9.8241—Il
U.
1.d52i—Lni U.
4. 4.16?+—us
3.3 j.
3. 4.1311—10
9.l07cLn3. 5.
5. 718+—uI 6.6501—li l.4631—sI 9.115+—nil .nnni 1. i3njc *uU 1.2Slni*oo
5. 5,~+—~i o.55.i+-o.,1 i.sOOo—si 3.4.?OL—C1 4.13+101,0 I. 149*. •so 1.4o3,.*,.,)
rn.9441—sl. 7.u? 14—ui 0.213+—ok 9. 3’,2.—niA A.o47n.*jo I • iou •n.’I i.2741*ni4
6.457+—ui 7. list—uI 3.3215—ui H.450*.—sI 1.0 *44
o.illt—u+
6.194+—IL
5.
1. jsJs—UL
?.4L5+—uA. 5.549+—nI 9.5311 —ni 4..n..30 1 • 133+ 045
7. 530*—UI
-no.333+—ui.
3. 4.1911—Is
—i.442*n—1
4.
Cm0
IAXIS
1. A7n~.030 i.~h)t*~u
2.4~0c—ni 3.Subl—Ji 4.574+030 A • £3 3,. *‘~o +.cO?L+oi
397n.—oi
39
F IT 4.it.7+—ui U.
u. 4.167+—UI
0.
4.407+—u 0.
n.~ 0.
...
—1.442+—,.4
FIT o. 0.
I. 4541—11
4.1041—UI
I.409E—i0
~.
o.27in.—nnni 7.iulc—ni1 3.4*6+—Il 9.1651—ui 1.4595*10 i.ulu+*uni i.1741*uAi.inOj++lu u.2~5i*uu l.2,7++ou
6.2440—41 ?.4101—uI 5.5494—UI 9.+0i*—~i i.441104U 1.1931*04
6. 3971—niL 7.530+—ui 8.6+21—Ui 9.715101 1.4931+00 1.2561+44
no.5101—ui 7.6431—Ui no.775+—UI 9.9081—uI i.1U41*U0 1.2171*44
6.6241—01
6.7371—UI
7.7561—UI 4.8811—Ui 1.UULLUUU 1.1154*40 1.2291*uo
7.4691—UI 9.0021—li 1.0131*60 1.1271*40 1.2401*00
1. 3.3,1+—is
0. 0.
1.454+—il
3.4?Si—iU
1.4091—14
0.
ni.
.,.,441—u,
4.1875—si
3.
U. 3.3’ist—Io
5.710+—Ui 5.450+—on. 7.963i—ni1 9.1151—UI 1.025u*44 1.i38c*uu 1.1511*00 39
5.031+—nil
5.,4+—s1 7.377+—SI o.2491—ol 1. 342+—si
o.u37s—Ui 7.LSr.nL—nnl 4. 32+1—ui 4.4o5i—ui
4.46l+—Ui
u. 4.ls7n.—~I
—1.o521—u1 U.
‘.
is 71—ni
0.
148
14015.
o.9#nj+—ll u.0561—i1 5.Ziot.—ui i.435+’sU I.i49co~J 1.253+003
I.347+.s0 1.1614*54 i.274~*1u
FIT U.
—no.3332—uI 0.
4.15714+ 0.
R, GruberetaL /ERATO NIFAC NTFAC
NIFAC OTFAC rITFAC
OTFAC 9TFAC nIFAC *ITFAI.
NTFAL OTFAC NIFAC ,*TFAC
NTFAC rIrFAC NTFAC
rITFAC N1F4C NTFAC $5181
SITFAC NTFAC
uTP&C
Nkl,38,NL.1RA,12
1,112,iLc N1,Z,rlO,NL,lrIl,121,142,It+
2
03 05
Nu02,Nk,NL.lnnI.ILi.,Il.s,12+ NOZ,Ni,NL,141.ILA,Irr(,IL4 N82,I1o,r,L,1%1,ILi,1nC2.12c
3
03
~
os
NRZ.N4,NL,I4l,1ui,1R4,IZ~
63
63
M82,HnI,.,t,icu.I2n.,IA2,l22 082,38,52,Iii.iZl.I*2,IiA NR2,9t,41,18l,11L,102,AZl NRZ,Iir(,n42,l1,L,ILL,Ii14,1u2 NI1l,N8,NL,111,LLI,Iil.112 NRZ,N*n,NZ.I,nI.IL1,1’2,1Zc
05
+3
o7 53 SI
4*3 uS
ru8Z,kA,N2,inrl,I2i,I+2.IZ.~
7’,
?u
N02.44,r12.14i.12i,154,121 N8L,MR,1,44.I1L,12n.,ilo2,122
7+ /3 i~
P,82,Nr,nMZ,Ioi,I21.1,n2,iZs NR2,N
is
5,NL,1k1,ILi,162,Ii.. NRZ,N*,91,1,*1,12i,I*4,122 482,N6,N2,1*l,I22,IIrd,I22 NRZ,N8,niL,III,IZn.,152,I14
NkL.r48,IIL,lli.1EA,102.112 N82,Nc,51,I,IA,I24,I’r2,12+ NR2,98,52,101.ILi,i42.1L4 nlR2,N0,n0L,irl,1L~.,t42,iL1
—.ou9111—01 —.173b61.uU —.144791,1,1 —.184481*10 —.171201+40 —.154941*00 —.640241—01 4
0:,
30
—.o73Lu+—~i .jsls.+—U+ 7+000 —.L3nn3,1*sn —.+no4n.41*su —.i/32 —.+449)t*u3 —.5+USoi,3i mu
47
0:, 0+
s:,
65 5, 53
33
+4
i
33 33 33 uJ
64 6’, 04 i4
2
33 33 33 ,3 33 ii 33 :,3 33 33 33
ii uI 30 2’I 4* 27
31
+0 +4
15 43 42
55 SS 05 sin
3 4 3
3c 32 32 ii 3i. 34 3~ 32 32
4 3
45 14 43 24
04
33 33
2
*1
1
01
4
7
III
03
33
0
i 1
i23 124 is,
+
12/
03
33 33 33 33 33
3
048
4) o5 03 55
2 2 2 4 4
73 —.354,11—04 0i.1u—ui —. +7os*.+,.uO —.l3:,721’sO —.102035*14 —. L4445c*30 —. 303932—.
3
33 33 33 33 33 .13 ii .13 uJ 33 33
1+u +41
—.
2
.12 ol 35 49 2+
1
no 7 6 5
4
£
4
3
4
3
i
4
+
2
1
1
1434+1—+i —.31353+—Ui .lT543lOUnn —.+40U5u*lu —.13125*054 —.+3*97+*niu —.235314—ui —.
369
i
03
17 25
stability code
—.16495+—Ui —.37903+—li —.140s01000 —.14)no2+*oO —.i79s9+oUs —.127o4+*s~ .O13?ns11c
3ZbilInnI —.254241—41 —.944604—UI —.162111*50 —.I5o:,3+*00 —.176421000 —.110441.00 .41U151—02
—. —.104330+01 —.18355t+uU —.lSSnnUEOuU —.A73431*Uu —.iUOi5E*Ul..
—.391731—01
—.468811—Ui
—.147051.v0 —.184541*00
—.
—.159221.00 —.169541*00 —.990651—01
Ilji3louu
—.185051*00 —.l63i6+*uU —.165221.0+ —.68481+—ui.
—.53930101.
—.ii9USE.00
—.185188*04 —.150831*011 —.160371*04 —.77u191—Ul
1+
FIT u.
3.54o+—u3
,.i57,.—UI
u.
ni. 4.1331—is
4.
0.
4.
0.
5.7101—4+ 6.450—si 7.9431—ui 9.1.151—Ui 1.o251*uu
5.83I~—S1 5.98.11—ui 5.04s1—ui l.2~nis—4L 1.0301*14 1.14*1+03 1.26,+.1.u
in.io4+—uL 7.077+—SI 0.lu+S1 9.34 c15 1 1.J47n.’ss 1.16150+0 +.274n.~1u
~.A47t—nii 7. Alit—ui no.32u1,u1 1. 4551—uI ~~1,591*o5 i.1721.so i.2uSt,00
0.Li++—Sn. 7.3431—UI 3.4,1,2—ui 3.10 ++~0n. i.ulU+*sJ 1.1331*00 1.4.7+0.5
s.204+—ul. 7.4061—Ui no.0.91—ui 9.5911—01 A.10i1*uO 1.li34.*sU
1.
4.333+si
4.Lb/+—ui
u.
0. 4,1331—lu
4.101+—niL u.
.n. 4.1572—si
3.1*44*—ni u.
.,.
4.,o71—0i
5.4311—si o.9831—ui. ,.u’1~i—s1 1.4235—ui i.53a1*nnl 1.149*055 1.2531*40
3.9441—uI
o. US7+~niI 7.1954—si 1.32~t—u1 9.4531—ink I.+:,9100s L.u?4s*Un.n 1.2451*oU
—1.8521—01
uc7L—uI
7u01—11
4.1841—01
2.4881—lU
0.
6.397+—Ui. 7. S3Uc—Uj 6.602+0I 9.795+—UI A.09310UU 1.2Uoi*UU
0.5101—UI 7.6431—Ui no.7754—01 9.9081—UI i.i041~0U n..2171*UU
6.6241—61 7.7561—UI 8.8891—01 1.0U21+uO i.i151*0J 1.2291.00
0.7371—UI 7.8691—01 9.0021—Ui I. 013E.UU 1.127E*uO i.24u100u
0. U.
—2.72,1—11
3.9811—11*
2.4891—l0
In.
u. 0.
—1. 7bU+—11
4.184*—UI
2.5361—10
u.
6.3971—LI 7.5341—UI *.6o21—01 9.7951.—li 1.0931*00 1.4151*1n4
0.5151—Ui ?.5431—uI ,.77~1—ul 9.9041—41 i.1041*uU 1.2171.00
6.6241—Ui 7.7561—~u1 6.8841—UI 1.0021*00 1.iI5l,UU I.2291*UU
6.7378—01 /.Ut,9101 9.0021—UI 1.0131*1*0 A.1271+uU 1.24U..Un3
3.
—2.
108
1..134tnosnn l.cSit’Us
AX1S~
39
FIT 4.157+—un. 0.
FIT —1.0541—Li 0.
3. 98
lUsU
0.
104
XIS~
5. 7lnot—Ui 6.OSuE—us l.9031—ui 9.1151—01 1.u231*uu 1.I38i*IU i. 2,11 *40 39
7.s71+—si. 1.4091—ui 1.344+—nil l.u47+*uu i.I5+1*niU 1. 274+ *sn3
silt—ui 7.3s3+—3i .3.43b1.n3n. 3.500+—nil i.07n12*i. L.1,~~0nU i • 2* l+*us +.
5.4*141—Ui 7.4101Ui no.349.—Ui ‘3.63+L—Ai 1.0011*00 1.1*51000
370
R. Gruber etaL
/ ERA TO stability code
JP~I
01*
OsIIT*
.252970~2 .s19i1943 .43235833 .323b356. .44859254 .4,224729 .36775902
.s777uo 7u
.1.1933887 .11*770044 .4110+534 .1)423914 .u159s7,.3
.01619807
• A~, A3inon.,5 .14400414 .+1924353 ..1397u97, .0010+445 • 1~0S35 ,o7 .~,U4ol+4 .0360145.. .s:,5js3:n+
.34244139
.32434 312
.258358~7 .07695517 .13904592 .01,741350 .34227044 .u77064ld .n.ni89.,177
0,14*5,0 .30448952 .1i5,314~ 1.43404410 +.*n7+Ijodn, i.,Ssl+.+, 3.nn3u3493s 2. ni9i43’,3:, i.9U3o4o’~5 1.24*/483+ .0in)~,,+ii .015344+0 • 4bL+531u .3o7,345no
.31235042 .49779144 .1 1743340
*70 .23403433 .2469*1431 .4,045223 ..1391*Dull 44nL97b00 .i49.41,n, • 1337*347 .1l1i5323 .2450)7+4 .44473+13
.s1734353 • oui’i SFn.o .01)09310 .5,3357104 .00290)40 .
04243171
3~S1 011*10 .35497573
.i30,..3u,. • 1447iu 15 • 440695 As
1.55341343 1.s313 *144 4.13734273
.94985319
.3117/3*2+
2.5)5L,,207
.036*444n.3
.27*14,40)
.05494+31 .351734951
.2743s..4, .14’llt,nns .5.17,/,43
3.3194*337 4.3, liiuic 3.337.7151 +...JuUt.11i
.1,54)8,1
.212o os41
..5 4nn~3u,~o
.n.15,000’0
...1o7’,191
1.13122.7+
.6)911744 .03099310 .u14u4s9S .uni4550o4 .u3i47 192
.i,79+SSo .+,+53/37 .i42ci,5, .43773+14 • 35+4443
.17j4+oz5
.0)104101
.5,4343.7
•351212~5
+0070333
.14234117
.133+473’, .++i70o/+ .44+79.04 .40171531
431110
.311351)7 .3334325’i .46 333i05.
.b018477i .35703045 .75903984
:,)94o
72
.35614/94 .21612533
.1243.301 .ni9I.0usI3 .v229YU*4
22 53io37 .09000+10 . :,r,64,’, *1 .3132,3/7 .llSuO4snl •
.4,nilnoin’1’, • 3177/37+ .3,11+0 + 3 • 33323014 .3)277018
.i)639*
:,u
*00
—.+n77,8751, 1303471
—.53393263
.s3403+uO 4j ,344 .s~i75)+4 .4110794’, .o..’n,n*.45i
—.75s4641*l —i.-1
i.24354300
—...
—ILn.Ln.n33334
—.33u,2:,2~ —.A5i*ils
.51)54070
—..,in.,4711n4 —.4iu74)/5
.un..?373o . ,~I 3,375
—.5471?i41 —.01.11,511
.s24n,73n—.04141240—. =
—.039307/3 —1.12795931 — 1.7,409908
—3.35069932 —5.307541*4 —4.o.b77~+I .03430347 1.11777126 .9432459k .05)24958 .8937s99~ .si9iu:,Sl .oI /o36o1 10307)70
.2
47+
4.01,52511’
.44+365*43 .3,533001
~i.lcs~4S4*.
.7421n339)
on’s
•
1,3+795
—.4344915, —.llIs+550 —.4)7,44c0
+1
s.+4o3741,5
4.
—.1,5333,57
—4.419753,’, —4.44271413 — 1.4*l+n+4j.~ .301*04277 .02533033 .7078)53+ .o2oio~,O • 4952093*1 .43036u21 .4)040794 .33+936u3
—2.01104157 3,147124 —4. 5007 i27 /..ni~711o9 —sO.31349308
4;’.j/.~~4
1. .oo4i/5/ 3.+3325i~4i 1.73301351 4.15,34,05 3.74 170590 4.42130 30,
8.34+.)no56 —5.iSsH13 —d.O4ouul3ln —3. .,so,94352 —2.us2o?1*53 —.77niuo93:, —.34.2,33*4 —. 17413)13
—.39886044 —1.40227441 —Z.j39159no6
.2305491, • 2s4.4737 • ci’,9143s .47043331 .,.20n.n.4u9 .u,720o4o .iuliUuq+ .,0
3/Si 019
*0*
—l.*1lo4n,HiO —4.24,37233 —3.1+7,3750 77
07*
~01540,
• 3205514
*00
l,uss’*/
.o
41314+0 .42450334 • 11444905 .s127 4S9i .u*423i .,9* • ~,21,418+ .ui2:,*7n,3 .45319134 .50494023
.01213205
.43541050 .2*3,suo7 .1372,0,9
.593,3174 • 74n,59+71 ./411+4,2
.004445+2
.335+0543
.70*3933+
.,u052343 .u405s4+u
.0o715211
.u0319775
.44+99+35
.s2i29602
.uu47nIsio—.
444
3.L
00*
ooun7inl
—.,.l,Su,h, ).977
191*31
—9.35045111
—43.bilOuiSO —io./6~7u764 —4.3+11*4,51 —1. 7004.23’u .00/0)+43 —.2354,’,91
—I. uo4ISo*49 —o.72030i11, —2.43*87055
—4.51174110 —7.j3so2lAA —4.49i33079 —. ~5
949434
2.44.35037 £ .o2204u 70
1.40344956
—.14303933 —.1,043)7)7
U5+3~36s —.0394114) —.03750+37
—.73131183
1.11438241 i.uU17ll144 .3,2831)5
.31411291 —.i737454u
R. Gruber et aL / ERA TO .3/SI *01070 .78370.43 .50340419 .64240694 .93524886 1.1955059b 1.01014193 .774u8i6o .498595,3
.3371.4+4 .42714344 .49054400 .,121.8943 .5I52’18u2 .5425434*3 .47490925 .43704463
1.19322413 l.u7inA*349*
.44481407
—14.4333957*)
3.443,?slb
4.l4,~oO3u4 3.23+24574 2.133,0947
.,.I9+077’, .1444443/ 759/ .nn33. .sno),l5l
—10.09991129 ~0./54)7832 —4.33881363
in.S)535044 i.lbioSmju .u5s3533, .04231542 .7o,~7o5i
.u24213o3 .0117501,4 .uu354:,,o .1u59?u?3 .ouu14~s9
.05071,11
.09453965
.017547u9
.40,o1733
.n,3.44355
.3321.377 • 47944*3, .5255025, .34533004 .43997745 .4,12.31,9 .40572,35
.141557*30
.40143275
.3370,063
• iS7sosY3 .2U5+,~~5 .15759,11, .12394313 .14253424 .u300304*, .o37niuS64
4IIX+1*4
—1.30129356 —.))5u5335 —.22471347 —.u9d15sU~ —.44072525 —.424u.408
—2.55153044 —3.49055+63
1.6560117? 1.74385126 ..32177047 1.31189787 1.1320+786 1.03615253
•03I~0jo/ .514o21.,ni—.
—.30450143 41103347
.9703911+ .94697030
.s1327+55
—.04979707
—.23622044
3751
1.0*1110913 .71946347 .7947444o 1.236nn1543 1.659059*3, i.529145uu 1.12502808 .73740044 .49059839 .3542Uc78 .28231034 .4394459? .41294985 .19819342 .19122520
—5.9755’8431 —0.91630780 —4.39U2999, .27487002
3.04979543
.23441101 .43919451
*00 —1.551.31163
2.44970044
.u72159o7 .06816438
*~)17T*
04*
3.973,7339 —4.73391330 —7.04934940
.1135+474 .175/4*4+ .47/3*7,3
.33357544 •3sn.,’9129+ .3149004, .4312908s .27i+7~o4
371
1.i
*70
.30871534 .20211,101 .141900)1 .10747535 .04534318
*1*
=
stability code
•
9*’4
1.47074414 .10354301 ~s4212z34 3.3383455’, 4.4,’,2544c 3.54317679 2.o0~1541 0.3742171) • ,71s 718+ .3/74)85+ .397.379*1’, .2+447923 .21I:,2nio0 • 173 + 7+nii .11227539
.u13741sn. • o45*3~L4 .160)4,73 .i9177023 .I343s993 •i35o2705 .o~no7005., .40204045 •uc4bOUis .014,4004 .u07U43oA .0543+724 .112,3430 • uu1b3399 .0ui365,4
.2=74?s04
.uL~55o54
*40 —4.07+87444 —5.8~7~2420 —9.~7/+o8d1 —15.89533192 —n.no.7715,135 —.u.1nnU39523 —3.93857077
—1.35/5in397 —.51340314 —.40.73023 —.0,027047
—.u3?o9913 —.51504552 —.337030.9 —.uo4sSSOo —.o)541595
*48 —2.31969036 —3.66089913
—~.5)420n373 —20.40016102 —11.97)12,73 —4.9550310o —.48092042 A.134u0229 1.15106791 .94050834 .~79lbI52 .011431,45 .51349331 .4,679041. .93139576 —.44,o*31i2
ERATO STABILITY CODE R. GRUBER, F. TROYON, D. BERGER
“,
L.C. BERNARD
~,
S. ROUSSET
~,
R. SCHREIBER
EURA TOM/Switzerland Fusion Association, Centre de Recherches en Physique des Plasmas, Ecole Polytechnique Fédérale de Lausanne, CH-1007 Lausanne, Switzerland
W. KERNER, W. SCHNEIDER Max-Planck-Institut für Plasmaphysik, Garching, EURA TOM-A ssociation, Fed. Rep. Germany
and K.V. ROBERTS Euratom/UKAEA Fusion Association, Cuiham Laboratory, Abingdon, Oxfordshire 0114 3DB, UK Received in final form 11 September 1980
PROGRAM SUMMARY Titleof program: ERATO
Keywords: plasma physics, ideal MHD, variational principle, finite elements, spectrum, instabilities, eigenfunctions
Catalogue number: ABVS Nature of the physical problem This computer code treats the stability of a Tokamak-like plasma described by the ideal linearized magnetohydrodynamic (MHD) equations [1]. The plasma is considered to be in an equilibrium state. After perturbation of such an equilibrium, the evolution of the normal modes can be calculated by linearizing the ideal MHD equations. Any mode (unstable or stable) of the spectrum can be examined.
Program obtainable from: CPC Program Library, Queen’s University of Belfast, Northern Ireland (see application form in this issue) Computer: CDC Cyber 170-720; Installation: Ecole Polytechnique Fédérale, Lausanne Operating system: NOS BE
Method of solution The variational form [2] of the 2D ideal MHD equations is treated by a finite hybrid element approach [3], which proves to be well suited to describe the features of the problem2Bx with sufficient accuracy 141. The eigenvalue problem Ax = w is solved by VEKIT (an inverse vector iteration), a subprogram of the block matrix library HYMNIABLOCK [5]. The code consists of 5 main programs (ERATO ito ERATO 5) linked together by disk files.
Programming language used: FORTRAN High speed store required: 35 000 words Number of bits in a word: 60 No. of magnetic tapes required: none Other peripherals used: line printer, disk files
Typical running time The running time is proportional to the number ofpoloidal flux surfaces Ni,, and to the cube of the number ofangular mtervals N~when N~>20. A typical case with 24 X 24 intervals in ‘I’ and x takes iS mm per eigenvalue on a CDC 6500.
No. of cards in combined program and test deck: 7901 Card punching code: EBCDIC
Present address: * IBM, Lausanne, Switzerland; ** General Atomic, San Diego, USA; *** Oak Ridge National Laboratory, USA; **** Zschokke SA Geneva. Switzerland. 323
Unusual features of the program The code uses subroutines from two program libraries, the
324
R. Gruber et al.
/ ERA TO stability code
utility routines of the CDC OLYMPUS package [6] and the eigenvalue package HYMNIABLOCK [5]. However, all the required subroutines are included in the ERATO package itself. ERATO is written in STANDARD FORTRAN [7], except for the use of the input facility NAMELIST which is available on most computers, and the CDC statement ENCODE which is used in ERATO 5. The plotting routines PLOTS, PLOT and SYMBOL used in ERATO S must be written to suit the local computer system and details oftheir arguments are given. ERATO also uses the CDC function subprogram DATE which picks up the actual date for plotting, and LOCF which gives the base address of an array. The programs ERATO 1—5 axe arranged so that they can readily be combined into a single program if sufficient storage is available. References [1] B.B. Kadomtsev, Reviews of plasma physics, vol. 2 (Consultants Bureau, New York, i966) p. 155. [2] W.A. Newcomb, Ann. Phys. 10 (1960) 232. [3] R. Gruber, J. Comput. Phys. 26 (1978) 379. [4] ~R.Gruber and F. Troyon, 3ème Colloque Intern. Méthodes de calcul scientifique et technique (Springer Verlag, Berlin, 1979); Lecture Notes Phys. 91(1979) 288. [5] R. Gruber, Comput. Phys. Commun. 20 (1980) 421. [6] M.H. Hughes, K.V. Roberts and G.G. Lister, Comput. Phys. Commun. 10 (1975) 167. [7] Standard FORTRAN Programming Manual, Computer Standard Series (National Computing Centre Ltd., Manchester, UK, 1970).
CONTENTS 1. Introduction 2. Equilibrium 2.1. General equations 2.2. Solovev equilibrium 2.3. Numerical equilibria 3. Stability problem 3.1. General axisymmetric problem 3.2. Choice of independent and dependent variables 3.3. Potential and kinetic energy in axisymmetric geometry 3.4. Vacuum energy in axisymmetric geometry 3.5. Additional n = 0 terms
4. Numerical approach 4.1. The marginal point of the continuous spectrum 4.2. The finite hybrid elements 4.3. Choice of the mesh 4.4. Spectral shift and its elimination 5. Organization and normalization 6. ERATO 1 6.1. Initialization 6.2. Choice of the mesh 6.3. Analytic treatment of the axis 6.4. Mapping 6.5. The normalized equilibrium quantities 6.6. Scaling law 6.7. Mercier’s criterion 6.8. Ballooning criterion 6.9. Definition of the betas 6.10. Surface fitting 7. ERATO 2 7.1. Treatment of the singularities 7.2. Discretisation 7.3. Vacuum energy matrix 7.4. Then = 0 terms 7.5. Resolution of the integrals 7.6. Definition of the shell 8. ERATO 3 8.1. Organization 8.2. Potential and kinetic energy 8.3. Symmetry conditions 8.4. Regularity and boundary conditions 9. ERATO 4 9.1. Eigenvalue problem 9.2. How to call VEKIT 10. ERATO 5 10.1. Diagnostics 10.2. Plotting 11. Test cases 11.1. Preset values (internal mode) 11.2. Test run output (kink mode) 12. Instructions to the user 12.1. Convergence properties 12.2. How to switch from equilibrium to ERATO 12.3. Choice of NAMELIST parameters 12.4. Control stream and data set 12.5. Dimensions 12.6. Space and time requirements
LONG WRITE-UP 1. Introduction Thermonuclear fusion research is profoundly concerned with the problem of macroscopic ideal MHD instabilities. When such instabilities occur they are generally very violent and global and lead to destruction of the confinement. The calculation of these instabilities has therefore become an important subject, and as a consequence of the great interest in Tokamaks, numerical stability codes have been written for axisymmetric toroidal geometry (2D).
R. Gruber etal.
/ ERA TO stability code
325
The first codes to be developed [1,2] solve the Newton equations of motion by finite differences. These codes work well for high pressure plasma but for weakly-growing modes the computer time becomes prohibitive and furthermore they are only able to find the most unstable mode. In the meantime two spectral codes (PEST from Princeton and ERATO) for a general axisymmetric equilibrium have been written and compared [3]. These two codes start from a variational form [4]. The existence of a continuous spectrum extending to the marginal point is very awkward for the numerical treatment. At this point, the mode is divergence-free and the B V operator vanishes. Weakly unstable modes, which are close to the marginal situation, are nearly incompressible and B V is very small. To overcome these problems, both groups use different finite elements for different components [5] to satisfy the incompressibility condition, and so prevent spectral pollution. The other condition (B V = 0) can be fulfilled by a Fourier analysis as implemented in PEST [6] or by a “finite hybrid element” approach [7] as used in ERATO. These hybrid elements have numerical advantages. They are piecewise constant in each mesh cell. The equilibrium quantities can also be chosen to be constant in each mesh cell; they only have to be given in the centre. The integration of the matrix elements is then trivial. The number of unknowns is smaller than for a regular finite element approach [8]. The convergence is generally from below, opposite to all Galerkin procedures. This enables us to well separate the unstable modes from the stable continuous spectrum and we never miss any unstable displacement. ERATO consists of 5 main programs (ERATO 1 to ERATO 5) linked together by disk ifies. For numerical equilibria, a sixth program, the equilibrium code, has to precede ERATO 1. ERATO 1 reads the data produced by the equilibrium code and maps them into ERATO form. It also prepares the analytic extended Solovev equilibrium [9] to be used as a test case. ERATO 2 calculates the vacuum contribution to the energy ~W. ERATO 3 builds up the matrices A(potential energy) and B(kinetic energy). The eigenvalue problem Ax = w2Bx is treated by ERATO 4 which makes uses of the subroutines of the eigenvalue package HYMNIABLOCK [10] included as part of the ERATO code. The eigenvector is analysed in ERATO 5. All 5 programs call subroutines of the CDC OLYMPUS package [16] which are also included. Subroutine maps are shown in figs. 11—16, calls to OLYMPUS subroutines being omitted. Tables 12—16 give lists of subroutines. The structure of ERATO 1 is similar to that of an OLYMPUS code but the order in which the initialization subroutines are called is not identical. The main programs ERATO 2—5 are dummies, and the entire code can readily be combined into a single program if sufficient storage is available. 2. Equilibrium 2.1. General equations In axisymmetric geometry, the magnetohydrodynamic equilibrium is given by VpjXB,
VXB=f,
17•B=0.
(2.1)
Writing the magnetic field as B=~ê~c,_±ê~,XV~,Ll,
(2.2)
we combine eqs. (2.1) and (2.2) into a second-order differential equation for /‘:
a /1 ~.,\ a2~ ~ -i--) = —j~r, where the toroidal current density4, is *
i/i
~--~—~-
dp
1 dT2
(2.3)
(2.4)
R. Gruber et al. / ERA TO stability code
326
A polar coordinate system (r, z, p) is used. The two functions p(i,Li) and T(i,ti) correspond to the plasma pressure and the flux of the poloidal plasma current, respectively. They are functions of ~tionly and can be arbitrarily specified. In order to prevent surface currents we require j~,(~i= t,115) = 0, where ~ is the flux at the plasma surface. The safety factor q (i~) is defined by dl.
(2.5)
The integration path follows a constant flux surface in the poloidal direction and B~is the poloidal magnetic field component given by 2B~= ViP2 = (~)2 + (~iP)2 (2.6) T
2.2. Solovev equilibrium A very simple polynomial solution has been found by Solovev [9] and is characterized by 4iP
— — —
2)(i~i iP~) 5(i a2E2R2 +E
~2
—
— —
~2 °
—
4iP~b2 a2R2E~i~,
27 ( ) .
leading to iP 2+b2 ~‘a2R4[5[r E2
2—R2)2 Z 2(r
.
28 ( . )
In expression (2.8) the quantity E denotes the ellipticity, R is the distance of the magnetic axis from the toroidal axis and a is the inverse aspect ratio. In this equilibrium_the toroidal current density j~,is flat. Note that on the symmetry plane z = 0 the plasma is in the range R.,~/f = 2a ~r
3. Stability problem 3.1. General axisymmetric problem ERATO solves the linearized ideal MHD equations [12] in variational form [4]: ~L=~(Wp+WV+Wk)=O.
(3.1)
W~,is the potential energy of the plasma: ~
fdr
VX(~XB)+(v~)(fX~)]2+yp(V .~)2 —2(v .~)2(jXv)(B~V)v},
(3.2)
R. Gruber et al. /ERATO stability code
327
where ~ is the displacement of the fluid elements and v the unit vector which points outwards normal to a p constant surface. W~,denotes the vacuum energy
wv=& fdr~v-
(iPViP)}.
=
(3.3)
0 is the scalar potential in the vacuum region. In the kinetic energy Wk
=
~-
f
(3.4)
drp(~-~)
p denotes the mass density of the plasma. One boundary condition demands that the component of the magnetic field perpendicular to each conductor c~vanishes =0. The continuity condition through the plasma surface (s) is written
(3.5)
Cj
[V X (~ X B)]~ (3.6) This general time and space formulation is reduced to the axisymmetric case by imposing homogeneity in the toroidal direction ~ (see fig. 1). This coordinate ~ as well as the time become ignorable variables after a Fourier analysis on the displacement vector =
~(iP~ x~p, t) = E~(iP,X)en~(~t).
(3.7)
Here n is the wave number in the toroidal direction and ~2 becomes the eigenvalue of problem (3.1). In that formulation, ~2 > 0 describes an oscillatory and stable solution, ~2 = 0 a marginal situation, and w2 <0 means that the equilibrium is unstable. z z
:___
a
=::::
-__
____
plasma surface
-2a AB
(Rest—i) a
Fig. 1. (s, x) coordinate system in plasma region (a) and (p, 0) in vacuum region (b).
328
R. Gruber et al. /ERATO stability code
3.2. Choice of independent and dependent variables Since in ideal MHD the fluid elements are frozen to the flux surfaces it is crucial to take 0 or any function of 0 as an independent variable. We found it convenient to choose as independent “radial” variable s =
where i,(i~denotes the toroidal magnetic flux at the plasma surface. The “poloidal angle” x is chosen in such a way that the Jacobian becomes J = qr2/T. In the infinite aspect ratio limit, J = J(0) only. This angle x has been used by the Princeton group [6] and proves to be convenient, since dX=q~ 11p1d1=j~-~1 ,
(3.8)
JB- Vfrr~/_+inqfEEqF(f)
has only a multiplicative equilibrium factor q = q(i,Li) in F. As dependent variables we choose the contravariant components 2,V 2siPs~x(rIViPI Y+t~3~~Xand Y=2s0 X=JIVi~’I~~/r 5~/Tr, —
parallel to the magnetic field. The quantity i3,~=
(3.9)
(axIa~)~ measures the non-orthogonality of the coordinate system.
3.3. Potentialand kinetic energy in axisymmetric geometry With the Fourier expansion (3.7) the variational form (3.1) becomes 2~K0. (3.10) 6L6Wp+~Wvw Introducing the independent variables (s, x) and the dependent variables (X, V, Y), the potential energy W~, becomes W~=1~ff~~~dx[aIIiI2 +b11 2 +c11 2 +d11 2 21 31 41
—
elI 5
(3.11)
2],
with
iax
~ ax av
11F(X)su——+inX, q8x
~
13 14
=
V ~ .
r1
2x aras
T 1 r
12=rVT—i=—+—, r as ax
a(Tx/q) + 2j~ir as lViPls ar2v q2 2
ax
~,
F(r Y),
+—~-
r
X,
I~=
(3.12)
and a
2q2i,lir4/J3IViPI2,
b
T2r2/20
2r2I2iP 5J,
c
IV’P1
4’y~/2iP 5J,
dr
8J,
4 4 ~-(lnlViPDn _.~~~(lnr)n]. (3.13) e_.±[j_ 12+ We have chosen this form in order to be able to treat the limiting case of a currentless plasma. This formulation has also analytic advantages since there is only one negative quadratic term which in a marginal situation must cancel the four positive quadratic terms. In a tokamak ordering (e = IViJ.i IT << 1), for example, the second term seems to dominate. Whenever ‘2 is not small, this term defmes the fast modes. For the other two classes of the spectrum this term must be small (12 ~ a for the slow waves and 12 ~ ~2 for an Alfvèn or a kink). Another advantage of (3.12) consists in the presence of Y only in one single, the fourth, term.
R.
Gruber eta!. / ERA TO stability code
329
The “kinetic energy” becomes
411 ~dX{fIXl2 +g(V+Y—f3~XI2+hlYI2},
(3.14)
Kr~ with
f= 2p~iTr2/qlV~il2, g = pIV~I2qr4/2T~
4Tq/2~ 5, h
pr
5.
(3.15)
In the plasma region we have to impose a regularity condition at the magnetic axis X(s0)0.
(3.16)
When we are only interested in internal modes we can put the wall right at the plasma surface. The vacuum contribution then tends to infinity and is replaced by the fixed boundary condition X(s1)0.
(3.17)
In addition we only look for symmetric displacements which are given by the condition (X, Y, V)(x) = (_X*, y* V*)(21r where
*
—
x),
(3.18)
denotes the complex conjugate.
3.4. Vacuum energy in axisymmetric geometry In the vacuum region the scalar potential 0 which appears in the vacuum energy (3.3) satisfies 20=0.
(3.19)
V By use of the vector identity =
V• (OViP)
—
0V20
(3.20)
and Gauss—Ostrogradski we obtain
(3.21)
Wv~ifO~ds+~~ fords,
where p and c represent the surface of the plasma and of the external conductors, respectively. For n = 0 and for a symmetric displacement, the compression of the field lines in the inner region of the torus leads to additional terms [131. On each conductor c• the boundary condition implies = 0. The continuity condition on the plasma surface (s) is
[VX(~XB)] Then W~becomes W~ —T ~
f
v=~ ~
r0~F(~ d~.
F(~.
(3.22)
(3.23) (3.24)
R. Gruber et a!.
330
/ ERA TO stability code
Including the boundary condition (3.5), we can use Green’s theorem 0(x)
=
~L f[o(x’ 2ir S
aG(x,x’) a~
G(x, x’)
—
~~11 ±~ ~ a~ ds’
—
fO(x’) aG(x,x’) ds’
2ir
(3.25)
ci
with G(x,x’) = 1/lx x’l, where x’ and x are the integration variable and a point, both on the boundary surfaces, respectively. In the numerical treatment it will be convenient to use polar coordinates (p, 0, p) to integrate (3.24). In this system: dz 2 i=~i~_, ds=lrd~,d0=rdl=JlViPldx, 3 l/dz a dr 3\ i/7dr’\2 —
~
+(_),
r= R +p cos0,
Ix —x’I
p sin0,
z
=
[(r
—
r’)2
+ (z
—
z’)2 + 4rr’ sin2 (p12)]112.
(3.26)
After having Fourier analyzed 0 and 3çb/&v on the surfaces 0(0, ‘p)l
5~= 0(O)I~,~e’~, 30(O,~ =N(0)e~’~’
(3.27) (3.28)
ir av and defined the p-integrals G~(O,0’) e1’~
2ir
,
G~n’(O,0’)=f ~
dp,
(3.29)
IX~X~’~
we obtain the two following coupled equations: 0~(0)=
~—fdO’ r’v’
0~(0)=±
f
do’ r’v’
V’G~’0~(0’)
V’G~ 5’Ø5(0’)
—-~—f dO’ r’v’- V’G~’0~(0’)
(3.30)
—~fdO’ r’v’
(3.31)
do’ G~’N(0’)
—
1f
dO’ G~’N(O’)
V’G~c’0c(O’)
—
with v
•
V
,~ dz’ 1aG~B’\
GAB’
—
~
\
a,~’J~,
—
dr’faG~B’~ az’ ~
ã~k
(3.32)
and A and B can both either be on the plasma surface (s) or on the conductor (c). The prime is set on the integration variable. In the present version of the code we restrict ourselves to one conducting wall only. 3.5. Additional n
=
0 terms
In the previous treatment of the vacuum contribution, we assumed that for the perturbed magnetic field in the vacuum region: jIiB. dl=~ViP- dl=0.
(3.33)
Condition (3.33) is always true for displacements with n 0, and for n = 0 and no m = 0 contribution. However, for n = 0 and m = 0, there are other terms which have to be added to W.,. due to a compression of the average
R. Gruber et al. / ERA TO stability code
331
magnetic field by the displacement. The perturbed magnetic field liB in the vacuum region can be written
6B=zi~+z2X~
(3.34)
,
where Z1 and Z2 are constants to be calculated, a fulfils A*ct=0
ac =0
a51,
(3.35)
and 0 denotes the scalar potential which has been treated before. The vacuum energy then becomes 2 ÷z~ (e~~X Va)2 W~= 4 J’dr
[z~(~)
(3.36)
Note that all the coupling terms disappear due to (3.33). In this section we calculate the first two terms of (3.36) for which a closed line integral in the vacuum does not vanish in all the cases of an axisymmetric mode. To compute the constants Z 1 and Z2 we form the integrals
f (dsp XA)
Ii
12 =
•~,
~
f(ds~XA)
(3.37)
where A is the vector potential corresponding to SB defined in (3.34). On the plasma surface d.c = rdp dI X ê~,. Using the boundary condition
(3.38)
(vXA)~=—~,•B, the integrals (3.37) become 11
=j~(dlXv)
B,
‘2
(3.39)
=j~,p(dlXv). (BX Va).
The integrals can also be converted into volume integrals over the vacuum region Ii=fdrV.(AX~!~L)’*’Zifdr3~,
I2fdrV.(AX~~)Z2J’dri-!,~_.
(3.40)
Here we have used the relation (3.33). By equating the expressions (3.39) and (3.40) we obtain the constantsZ1 and z2 in terms of ~ on the surface. By substitution into the expression (3.36), the final vacuum energy becomes 2/fdr-~ + [j~@1xv). (BX Va)]/ fdr~!~_}. (3.41) W~= 4([1~(d!x v) B]
Knowing that —4n~’Gk~’ is the Green’s function related to the operator
1=~
~-~-~-
_~-_J
v’. V’(—4rr’G~~’) —
~-
v’~V’al’(—4n-’G~ 5’)+
we can write Green’s theorem
~
~J— I
V’al’(—4rr’G~’)
(3.42)
~—v’. V’al’(—4,r’G~’)
(3.43)
‘•
on the plasma surface and 0
=
i-_f
~-
v’- V’(—4,r’G~~’)
on the conducting wall.
—
f
~
v’. V’~’(—4n~’G~’) ~
f
332
R. Gruber et al.
/ ERA TO stability code
In (3.30) and (3.31) we knew ViP and had to calculate
0. Here, we know a and have to calculate Va.
4. Numerical approach 4.1. The marginal point of the continuous spectrum In the 1D stability code THALIA [14] we had to choose our basis functions such that we were able to reach the marginal situation (w2 = 0). This spectral point was characterized by the condition V•~=0. (4.1) Spectral pollution [5] then disappeared and the stable part of the spectrum proved to be well separated from the unstable one. We adopt the same idea to solve the 2D axisymmetric stability problem and start to calculate the conditions defining the marginal situation in toroidal geometry. They are given by the three identities [15]:
(4.2)
F(V)=-~-+inV’0,
q ax
(4.3)
2
r
3s
r
—
3x
+~j- F(Yr2)~’0.
(4.4)
r
Eq. (4.2) suggests a Fourier analysis in x~[6], eq. (4.3) imposes the choice of basis functions in s and eq. (4.4) simply defines Y. Another possibility to fulfil eqs. (4.2)—(4.4) is the use of finite hybrid elements [7]. 4.2. The finite hybrid elements These involve extending the variational formulation and making a finite element approach in s and of solving ~L(X, V, Y)=0,
x~Instead (4.5)
we vary over 7 variables and solve 6L(~-~-_--, ~(2)~~T~
~
j~(2)
~
y(2)) =
0,
(4.6)
restricted by the evident identities ~1)
=
=
~
ji(l) =
v~2~ y(l)
=
y(2)
(4.7)
The original problem is reobtained when the identities, eqs. (4.7), are substituted back into eq. (4.6). In our approach, however, we vary the Lagrangian over all 7 variables ~ ~ ~(3) ji(1) ~(2) Y(’) and y(2), and do not impose the identities in eqs. (4.7) everywhere, but only on specific points in each mesh cell. To prepare the discretisation of the problem, we rewrite the constraints (4.7) in the equivalent forms urn
—~-
~
f(x(2)
f(v(2)
—
X~’~) dr = 0,
V~E ~l,
—
v0)) dr = 0,
V~E ~
lirn
-~--
~
f(x12)
f(y2
—
x~3~)dr = 0,
Vi~E ~2,
—
Y(1)) dr = 0,
V~E ~l,
(4.8)
R. Gruber et al.
/ ERA TO stability code
333
where ~l denotes the plasma domain. The conditions (4.2) to (4.4) then become:
av~1~
(4.9) =
+
0, 3lnr2
~)
~2)
(4.10) + ~2)(inq
(4.11)
+3~2)+~~~=o.
4.3. Choice of the mesh In order to discretise the problem we cover the domain ~l with a rectangular mesh. In the radial direction we choose N~intervals (the radial mesh point counter i varies from 1 to N 1 + 1) and in the poloidal direction we choose I / ~ N,, + 1. All the mesh cells can have an arbitrary size. In order to be able to impose the symmetry condition (3.18) we add two half mesh cells in the plane ir ~ x ~ 2ir (see fig. 2). Then we can confine the solution to the half plane 0 ~ x ~ ir only. We now satisfy the constraints (4.8) by identifying ~ as the size of a mesh cell L~,j,such that 2)—X~3~)dsdx0,Vi,/, I (X~ —~‘~)dsd~0, Vi,j, (X( ii ~ij ii ‘~
±
f f
~—
(V(~~ —
V’~)
ds dx = 0,
Vi, j,
_L
ii ~ij
f f (y(2)
—
Y(1)) d.c
dx = 0,
Vi, j.
(4.12)
‘I
These conditions do not uniquely determine the elements, so that we still can choose the order. We make the simplest choice possible by requiring that each argument in eq. (4.6) be constant on each mesh cell. Calling d 1 the basis function1’for a11 the basisthefunction for can ~ beg1+1~2,~ the and~hi+l/2J, thebasis basisfunction functionfor for~v<2) f~ andthe y(2), expansion written: basis function for V~’~ and Y~ NS+1 Nx+l
NS+1 Nx+l = ~
~
2~ = I~
X 11d11,
i=1
j=1
x~
i=1
NSfl Nx+l
~
j=1
2)
1=1 j1
(V) g~+112,1, Y 1+1/2,1
11j~1,
x~
1~11,
NSN+1
(V)(l) =
x
3~ =
x,
~
~
= N~N+1
(V)( i=i 1=1
(1
i=1
Y
1+1/2,1
j=1
h
1+1~2,1,
(4.13)
where the positions of the nodal values in a mesh ate all shown in figs 3 and the shapes of the corresponding basis functions are shown in fig. 4. It is easily seen that the expansion (4.13) together with these basis functions fulfils 2) reduce to the requirements (4.12). With this expansion the 7 unknowns ~ ~ ~(3) 1/(1) 1/(2) y(l) and y( the 3 unknowns per mesh cell X. 1, V1÷112,1and Y,~112,1. Note that d11(x) and ~ are functions linear in x and piecewise constant ins. The x derivatives of such elements are piecewise constant. A similar roof function linear in s and constant in x is used for f~1(s)such that its s-derivative is piecewise constant. The bases for f~1and h1+1~2,~ are piecewise constant. The basis functions for the radial component X are non-zero over 4 mesh cells, whereas those for V and Y are non-zero over 2 mesh cells only. The conditions (4.9)—(4.1 I) then can be fulfilled when they have to be fulfilled. We can now easily write each term in the Lagrangian (4.6) at the mesh centre:
R. Gruber eta!.
334
I ERA TO stability code
Ix
x
j 6
— .~_.
x~,
——
1,J.1
VL,lh 4
—
3
—
j.l
Nx5 N5=4
2
s
0—1
2
3
45
VL1/2j
Fig. 3. Position of the nodal values in a mesh cell for the “finite hybrid elements”. The physical values are situated in the centre of the mesh cell (o).
-
Fig. 2. Discretization of the domain l~.
ax~+112,~+1,2 x~+1,1+~ +x1,1+1 — ~ 2(Xj+i —xi)
X~,i.1
—
ax
x~2~
_x~~÷x
1,1+1 +X~+1,1+X1+1,1+5 4
—
_____________
—
v
v
Ai+l,j+l
Af~1,/— A1,j+l
2(s,+i av~+~/2,J+l,2 = V~+112,1+1— —
__________=
Yi+l/2,f+l
—
—
V1+l/2,/ x~
1/(2)
=
1+1/2,/+1/2
‘
2)
y(
~‘i+l/2,/
—
Au
s~)
=
V1+j/2,/+l 2 Y
+ V1+112,J
1+i/2,/ +
Yi+l/2,i+l
(4 14)
i+1/2,j+1/2 2 These relations can be considered as centred finite differences applied to a variational problem. Since in the Lagrangian, eq. (4.6), all the functions and their derivatives are piecewise constant, the equilibrium is also taken to be piecewise constant. This is done by replacing all coefficients coming from the equilibrium by their values at the centre of each cell. With (4.14), the boundary conditions, (3.16), (3.17), take the simple form —
x~,I=0, xJ%J1f1,J=o,
‘
l~/~N~+1,
(4.15)
and the symmetry conditions (3.18) in x yield X~’—X~, X~1=X~2, 1~i~N1+l,
~ = ~ ~ 1 = ~+l/2,2 V~÷1121 =—V~+112,2, Y~+l/2,l= Y~+1~,2
1
~ i ~N5.
(4.16)
We can now verify that all the three conditions (4.9)—(4.1 1), being piecewise constant, can indeed vanish over each mesh cell.
R. Gruber et al.
/ ERA TO stability code
7
j+1
~
335
d
e
j
.1~
~
1 ~
x
j+1 ~
~
+
g
~
1/2
Fig. 4. Shape of the basis functions d 2~, of f~ 3), ofg,~ 11(x)0)for and they(l) expansion and ofh~+jI of XW, of e11 for the expansion of 2)ax( nd y(2)~ 1for the expansion of x( 1j2,fom the expansion of V 2,jfor the expansion of v~
336
R. Gruber eta!. / ERA TO stability code
Table 1 The 5 ERATO main programs Program
Contents
ERATO 1
initializes, chooses equilibrium, interface equilibrium, stability, Mercier’s criterion, ballooning criterion
ERATO 2
calculates vacuum energy contribution
ERATO 3
calculates potential (A) and kinetic energy (B) matrices
ERATO 4 ERATO 5
solves eigenvalue problem Ax = o.~’ printing, plotting and display of the eigenvalue, the eigenfunction and other diagnostics
2Bx
4.4. Spectral shift and its elimination Compared with the usual finite element method [8] the additional freedom introduced by the enlargement of the functional space leads to a lower energy level. This helps to separate better the weakly growing unstable modes from the stable continuous region. However, for the Solovev equilibrium we observe a spectral shift [15] due to the numerical approximation of eq. (4.9). This spectral shift can mostly be absorbed by a small modification of nq. For numerical purposes it is convenient to impose the shift correction on the wave number n. This means that instead of calculating with an integer n value we introduce such that ~‘
(4.17)
i~n/{l —a(nq/N~)2—I3(nq/N~)4]. For a Solovev equilibrium we found a = 0.822 and 13 become smoother.
=
0.142. The convergence properties for fixed nq then
Table 2 Disk units Unit
Disk file name
Created
4(MEQ) 5 (NIN) 6(NPRINT) 7(NDA)
TAPE 4 INPUT OUTPUT TA
ERATO ERATO ERATO ERATO
8(NSAVE)
1 1 1—5 3
Used
Contents
ERATO 3,5 ERATO 1 ERATO 1—5 ERATO 4
mapped equilibrium read NAMELIST from card line printer unit matrix A NAMELIST on disk matrix B matrix LDLT numerical equilibrium scratch unit shifted A eigenfunction plotting quantities scratch unit surface equilibrium vacuum matrix
OUT
ERATO 1
ERATO 2—S
9(NDB) 10(NDLT) 11 (NDORY) 14~NSCRTC) 14(NDS)
TB TAPE 10 DORY TAPE 14 TAPE 14
16(NDES) 1 7(NDMAP) 17(NVAC)
TAPE 16 TAPE 17
ERATO 3 ERATO 4 equilibrium ERATO 3 ERATO 4 ERATO 4 ERATO 1 ERATO 1 ERATO 1 ERATO 2
ERATO 4 ERATO 4 ERATO 1 ERATO 3 ERATO 4 ERATO 5 ERATO 5 ERATO 1 ERATO 2 ERATO 3
R. Gruber et a!. /ERATO stability code Table 3 Common blocks No.
Name
Used by
Contents
C2.1 C2.2 C2.3 C2.4 C2.6 C3.1 C3 .2 C3.3 C4.l C5.1 C9.1 C9.2 C9.3 C9.5 N
COMPHY COMEQU COMEIG COMAUX COMVEC COMESH
ERATO 1—5 ERATO 1,3,5 ERATO 1,4 ERATO 1,3 ERATO 4 ERATO 1,2,3,5 ERATO 1—S ERATO 4 ERATO 1—5 ERATO 1—5 ERATO 1 ERATO 2 ERATO 3,4 ERATO 5 ERATO 1—5
physical quantities equilibrium eigenvalue constants auxiliary quantities eigenvectors mesh numerical quantities eigenvalue constants control values disk units mapping vacuum matrices A and B plotting quantities namelist
COMNUM COMIVI COMCON COMOUT COMMAP COMVID COMMTR COMPLO NEWRUN
5. Organization and normalization The ERATO code consists of five programs ERATO 1—5 which are explained in table 1. In addition, three external packages, namely: the numerical equilibrium (actually we use the ORNL code [11]), giving all the equilibrium and geometrical quantities to the stability calculations, some utility routines from the OLYMPUS library [16], the eigenvalue solver HYMNIABLOCK [10], are employed. All these codes are linked together by means of disk files which are described in table 2. All the disk input—output is performed in the subroutines IODSK 1—5. Table 3 contains a list of all the COMMON blocks used in OLYMPUS, HYMNIABLOCK and ERATO 1—5. —
—
—
Table 4 Normalization Position
rph = RTE
Displacement
tph = RtE
Fields
Bph = B0BE
Mass density Pressure Frequency
Pph = POPE 2o/12o)PE ‘~ph= (B ~ = (B~/(poR2~o))w~
Current
‘ph = (BORI,2o)JE
Notes: R is the real distance of the magnetic axis from the symmetry axis, B 0 denotes the toroidal magnetic field at the magnetic axis, P0 is the mass density at the magnetic axis, ‘E is the total current which is calculated and printed out in ERATO I (see section 6.9), eqs. (6.21).
337
338
R. Gruberetal. /ERATO stability code
At the end of section 12 we give a subroutine map of each main program ERATO 1—5 and of the complete code. Table 4 contains the normalization of the equilibrium quantities, the eigenvalue w2 and the eigenvector ~. Table 5 Namelist NEWRUN Quantity
Type
Default value
ALARG ALO ANGLE(16) ARROW ASPCT B2R2 CPSRF CST ELLIPT EPSCON EPSMAC GAMMA QIAXE
R R RA R R R R R R R R R R
0.0 —1.5 0.0 1.0 1/3 0
REXT SCALE WIDTH WNL WNTORE ITEST LENGTH MEQ NAN NCHI NCOLMN NDA NDB NDES NDLT NDS NFIG NIN NITMAX NMESH NPRNT NPSI NR NSAVE NSCRTC NSHIFT NSUR NTCASE NTURN NVAC NZ NLTORE NLY
R R R R R I I 1 I I I I I I I I 1 I I I I 1 I I I I I I I I 1 L L
4 10~ i0’2 5/3 3.3333 0.0 1.0 1.0E5 —2.0 1 4(Nx 4 ~JN5 15 8(Nx 7 9 16 10 14 0 5 10 0 6 14 65 8 14 0 15 2 10 17 33 T T
+
l)(8Nx
+
1)
+
9)
paper size for plotting eigenvalue guess angles at which displacement plots are performed length of arrow inverse aspect ratio a b2R2 for analytic equilibrium set in ERATO 1 set in ERATO 1 ellipticity squared E2 accuracy of eigenvector machine accuracy specific heat ratio inverse q-value on axis 1/q 0 defines conducting wall scaling factor parameter for non-equidistant x mesh dummy parameter toroidal wave number n CDC plotting parameter defined in AUXVAL channel number for equilibrium defines region to fit on number of poloidal intervalsNx defined in AUXVAL channel number for A channel number for B channel number for plotting channel number for decomposed A scratch unit number of displacement plots input channel number maximum number of iterations choice of mesh output channel number number of radial intervals N~ number of r intervals in equilibrium channel number for NEWRUN scratch unit shift in q plasma surface for plot selects the case number of turns for ballooning criterion channel number for vacuum number of Z intervals in equilibrium T = toroidal geometry T = compressible case,F= incompressible case
R. Gruber et al.
/ ERA TO stability code
339
6. ERATO I
6.1. Initialization ERATO 1 first reads the 4 label input cards. In CLEAR all the COMMON variables should be put to zero. PRESET then sets the first test case later discussed in section 11. The NAMELIST input (see table 5) is read in DATA and in AUXVAL we calculate some matrix sizes used in later main programs. Then FORMEQ selects all the cases defined by the NAMELIST quantity NTCASE, described in table 6. For NTCASE = 3 or 4 numerical equilibrium quantities must be transferred via disk. This is foreseen in IODSK1 fork = 6 (table 7) which reads the number of mesh points in rand z, the equidistant rand z mesh 0(r, z), and, depending on the numerical equilibrium code, constants to defme the arbitrary functions p(~) and T(0). For our equilibrium code this statement is
READ(NDORY) NR,NZ,((CPSEQ(JR,JZ),JR 1 ,NR),JZ=1 ,~NZ),(EQR(JR),JR1,NR), (EQZ(JZ),JZ= 1 ,NZ),CPSRF ,(APP(J),J=1,1 0),(ATTP(J),J=1 ,10)
(6.1)
The way one reads NR, NZ, CPSEQ, EQR, EQZ, CPSRF, APP and ATTP is free; it is only necessary to be consistent with what has been written on disk in the equilibrium code. The two arbitrary functions p(~) and T(0) have to be defined in TANDP and must agree with the functions chosen in the equilibrium code. They have also to be normalised such that they correspond to eq. (2.3). In our equilibrium code the chosen functions are
TT’(O) =
=~
~
ATTP(J) ~‘.
(6.2)
The quantity Po is defined by the condition p(~ = 0~) = 0. The normalization with respect to R and T(0 is performed in subprogram NOREPT. The plasma surface is defined by passing the value i~ion the surface, i.e. CPSRF(tL’5) through (6.1).
0)
6.2. Choice of the mesh The existence of singular nq = m (m integer) surfaces which are concentrated at the outer plasma region requires a careful choice of the mesh in the s-direction. This is done whenever NMESH * 0. First we go through the whole mapping procedure with an equidistant s-mesh in order to count the number NSING of singular surfaces for which nq = m.
Table 7 IODSKI(k) k Table 6 NTCASE ______
21
rewinds all disk files writes NEWRUN on NSAVE
NTCASE ______________________________________________________
3 4
writes EQ on MEQ writes plot quantities on NDES
1 2
analytic Solovev equilibrium analytic equilibrium through interface
3
numerical equilibrium through interface fitting of the plasma surface
5 6 7
writes vacuum quantities on NVAC reads equilibrium from NDORY not used writes surface and labels for plot
4
____________________
8
340
R. Gruber eta!.
/ ERA TO stability code
In a second step, the distribution of the mesh points in the s-direction is controlled through a monotonic weight function W(s) satisfying W(0) = 0 and W(1) = 1. The locations of the mesh points s~,i = 1, N5 + 1, are ...,
given by s~= W’((i
—
1)/N5), ‘c/j,
(6.3)
1 (x) is the inverse weight function.
where s= W x Designating
1, / = 1, NSING to be the locations of the NSING singular surfaces lying within the plasma, the density of mesh points is chosen to be 3-i-NSING NSING 1 2+~ (6.4) ds3(1+NSING)~ 1=1 (s—x1) For the width L~we use = [ndq/ds];~.The normalization constant a is determined by W(l) = 1. The weight function then becomes
dW
3+NSING W(s) = 3(1 + NSING)5
NSING +
a
j=1
arctg
s—x~ +
arctg
x ~—
,
(6.5)
where 2NSING a=
1NSING ____________
3(NSING+1)i
‘
1=1
arct g
1
—~.
+
arct _.1~ g~
66
When NSING = 0 the mesh remains equidistant. The constant term in the density (6.4) prevents too strong a concentration of points around the singular surfaces. The formula chosen does not pretend to be the best one. We have sometimes multiplied the ~j’s with a factor ~ to put more points around the singular surfaces. In the angular direction x we have chosen an equidistant mesh. It has appeared that it is sometimes advantageous to concentrate the x’~near the outer plasma region, i.e. around x = 0. This is specifically true when a ballooning type of mode is to be represented. The user is free to modify the subroutine MESH. 6.3. Analytic treatment
of
the axis
2z2) around the magnetic axis we perform a 10 parameter fit To obtain an adequate knowledge of ~/i(r ~L’ —~f~ +f 2 +f 4 +f 6 +f 2 +f 4 +fjoz6 +f 2z2 +f 4z2 +f 2z4 (6.7) 3r 5r 7r 8z 9z 12r 14r 16r over (NANR—NANL)*NANZ equilibrium points. This region around the magnetic axis is defmed by NAN which denotes the number of s surfaces which should lie in the “analytic region”. The default value is NAN = ~ The magnetic axis R = RAXIS is found by solving ,
aiP/ar(z=0)=0=f
2+3f 4. 7r
(6.8)
3+2f5r 6.4. Mapping
We assume that the numerical code furnishes ~ (r, z) where the r and the z mesh are equidistant. For the stabi-
lity code we have to invert the information, i.e. for a given tji we have to find the corresponding r and z coordi-
nates. This is performed in subroutine INTFAC. We start our treatment at the z = 0 axis in the outer region of the plasma and first fmd r(iP, x = 0) by a linear interpolation between the two neighbouring radial mesh points. By means of the subprogram ADVANC we
R. Gruber etal.
/ ERA TO stability code
341
~
r(s,X=lr)
R
r(s.X~o)
Fig. 5. The 6 points which define the polynomial (6.9) in the mesh cell (r
1
advance in the positive x direction giving for each mesh cell the two ingoing points (nos. 2 and 4 in fig. 5) between which ~,1ilies, and searching the two outgoing points (nos. 4 and 5 in fig. 5). The spatial of the 1 = 2 anddependence z2 constant line is, for each mesh cell, approximated by a third order polynomial in r ~i(r2, z2)”a 2 +a 4 +a 2 +a 4 +a 2z2. (6.9) 1 +a2r 3r 4z 5z 6r To determine these 6 parameters a 1to a6, this polynomial goes through all the 4 edge points (nos. 2 to 5 in fig. 5) of the mesh cell. The two additional points are chosen to lie inside the surface and so that there are at least 3 distinct rand z coordinates. When we cross the z = 0 plane again, the angle x (3.8) must be equal to ir which defines the q-value of this =~/‘ surface (2.5). The ~,1iderivative of q is given by the periodicity condition on the non-ortho13~(j3,~(O) I3~(7r) = 0) gonality 1 .-i’r i /~ç~~2 dj3~= dx ~j-~- i~i~ (6.10) -~~-+
1.1 ~-
—
~—~-~----
— —
which leads to
dqq t’d
j
/~,r
X
1 dT +
1
Iv~I2~
611 .
(.
This simple procedure is very accurate as long as we are not too close to the magnetic axis which is achieved by the 10 parameter fit (6.7). 6.5. The normalized equilibrium quantities In the subprograms EQANAL (analytic treatment) and in NUMEQU (numerical mapping), we fill up the array EQ(20,NCHI) which contains all information about the equilibrium for an s constant surface. The 20 quantities in EQ per mesh cell are explained in table 8. All the EQ(20,Nx) arrays are stored on disk file MEQ = 4. There are N 5 such arrays on that file. In all the quantities defined in table 8 we assume that the flux function Ton the axis is normalized to 1 and that the distance the toroidal axismeasured r is normalized 2 thatfrom is obtained is then in unitstoofthe distance ofthe magnetic axis R (see table 4). The eigenvalue w w~B2 2, (6.12) 1]popR where BT is the toroidal magnetic field at r R.
342
R. Gruber eta!.
/ ERA TO stability code
Table 8
EQ(k,i),j1,Nx
-_______________________________________________________
k 1
s-coordinate of lower left hand edge
2
x-coordinate of lower left hand edge
3
s-coordinate of upper right hand edge
4
5,6
k-coordinate of upper right hand edge (s, x) coordinate
7
mass density p
8
plasma pressure ‘yp/qo~l
9
flux
10
free
11
q/qo
12 13
poloidal magnetic field ~pr non-orthogonality i3~= (axIas)~
14
distance from toroidal axis r2
15
(a lnr2/as)~ (a In r2/ax)
16
5
function T
2fqolV~ 2
5
2
17
~
18
K=—I
slVVd
T ds \q q
2~r ~2
/
dp a
a
___~_~+_~!.
—(In IV~I)~———(Inr)
qoLlV~PI r a~p
dXL’
2/J=T/q 19 20
a~
0
rfree
6.6. Scaling law By means of a scaling factor a (in the NAMELIST input SCALES), it is possible to transform one equilibrium into another one [23]. Calculating one numerical equilibrium makes it possible to perform a whole stability study in q 0 as a parameter. At the same time we set to 1 the toroidal magnetic field at the magnetic axis T0 = 1. The scaling law is given by 2p, ~—~-al/i, ~~-+aB~, ~-~‘a aT~, D2 T~+ a(T2 Tg), ~ -~q2~ (6.13)
4-+a/~,
~
-
-~
6. 7. Mercier’s criterion
Mercier’s criterion [26] can be written in the form M=~((
1)2
2(l+q)
_P’T~(Ci+C l+q
2)÷P1Cs(C4+T2Ci)+p12{~(~_C1C3)_C4(C1+2C2÷C3)]}~0. (6.14)
R. Gruberetal. /ERATO stability code
For a given ~-surface, the integrals C1, C2, C3, C4 and C5, with 2/T, J’qr
343
(6.15)
become
q
~
d~
Ci~J
—
C
~
q
~
(i2
—
R2)2 d~
2~J
2
22
(r —R) dx
C42ir~,
C5
rf
(6.16)
~
This form has been incorporated in subroutine MERCIR which also uses INTGRA for the ~-integration and QUADRA for the parabolic fit. As a result we getM~0 for a stable situation, while M <0 predicts an unstable localized mode. 6. & Ballooning criterion On each constants surface we test in subroutine BALOON if the modified potential energy [17,21] = ~-fdx-~-
(A
I~f2~+BXoI+ CIX0I2J
(6.17)
is positive definite (ballooning stable) or not (ballooning unstable). The quantitiesA, B and C are given by [see (3.13) for the definition of a, b, c and e] q’(x XoY~ q / abc ( 2jVir T’\2 C=(b+c4~\sIVlpI2+~~:;/ —e.
A
=
a
+ —~--I’13 +
—
B=
‘
bc
(13
(b+c)A\”~
+ ~‘(x
—
q
T2/’ Xo)~(2/iPr + T’\ /~~V~p(
(6.18)
These expressions have been obtained after expanding X and V for large n as
~
(x 0
v= (v0 ~
~
(6.19)
As periodicity condition we impose
Xo(~= —lir)
=
(6.20)
Xo(~= iir) = 0,
choosing ito be a large number (up to 100). The printed array BALLCRIT contains the number of negative eigenvalues of the 5 W~,matrix where the values 0 mean stability and the values >0 instability. Close to the magnetic axis 1 and inprecise information can be produced. However, we know that around the axis MERCIER and BALLCRIT are equal. -~
6.9. Definition of the betas Let us define the integrals (in BETAST) 2, I~ di~d~p, I~2=JfJd~d~p
~ff~
1B2
fJ”1 d~d~,
344
R. Gruber et a!.
zv-]fJdiidx,
/ ERA TO stability code
IE=f~f~id;l1dx/~,,lRff~d~l1dxp.
The usual beta value (/3), the poloidal (j3~)and the Princeton beta (/3*) are given by (3 = 21p/1B2, i~1= 81rIR/I~, (3* = 2~/’l~i~/IB2.
(6.21) (6.22)
6.10. Surface fitting The calculation of a numerical equilibrium takes a lot of computing time whenever discrete conductors are prescribed. In this case, the resolution of the equilibrium is not good enough to give accurate information about stability. To eliminate this problem, we pass twice through equilibrium codes. We follow the steps (I) Calculate the numerical equilibrium prescribing positions of the wires and the currents flowing in them (ORNL code). (II) Fit the position of the plasma surface (NTCASE = 4). (III) Solve the numerical equilibrium prescribing the plasma surface as found in (II). The resolution of such a code is much better and, in most cases, enough for stability studies. In this section we described the part of ERATO 1 which performs step (II). In this step we keep the polynomial factors of eq. (6.9) for all (r, z) mesh cells through which the surface curve s = 1 goes. For practical purposes we propose to use an efficient free boundary equilibrium code including a “Buneman Solver”. Such codes enable us to resolve the equilibrium problem with sufficient precision (up to an (Nr, N~) mesh of 512 X 256) [24]. An ideal way to resolve the fixed boundary problem would be a finite element approach in which the nodal points are iteratively put in the same ~(i = constant surfaces that are required by ERATO. Such a code already exists (SELENE [25]) and will shortly be combined with ERATO. 7.ERATO2 7.1. Treatment of the singularities In this main program we calculate the vacuum energy contribution expressed in terms of the normal displacement X on the plasma surface. In eq. (3.21) we replace 0 and aO/ap by X using eqs. (3.23), (3.28), (3.30) and
(3.31).
Defining the modified elliptic integral K 0 to be
1/2 [cos~~~n~]1/2
K0~~) =
(7.1)
and
2 + (z~’ zA)2 ]/[(r’a’ + rA)2 + (Z’~’ ZA)2], [(r’B’ rA) the real and the imaginary part of the integral (3.29) become
(7.2)
G~~’ =4K
0/[(r’~’+rA)2 + (z~’ ZA)]. 0(~)(—1) The integral K 0 as well as its derivative dK0/dr! exhibit singularities, since for small n
(7.3)
Kn(~)=(_l)0~1(~)+...,
(7.4)
—
—
=
—
= (_1)0+1
—
±~+ ....
(7.5)
R. Gruber et al. / ERA TO stability code
345
The logarithmic singularities (7.4) in the source terms of the two coupled equations (3.30), (3.31) will later be treated analytically. It turns out, however, that the poles in the integrals containing u’V’G~~’ must be studied carefully and, numerically, be treated in a consistent way. Only when this is done properly are we able to move the wall continuously towards the plasma surface [22]. In our case, we treat these singularities by adding and subtracting the same types of integrable singularities, namely d0’r’v’~V’G~’= 0,
2ir
~ ~-
~—fdO’r’v’
dO’r’v’’ V’G~°5’ dO’ r’v’ V’G~~’ = —2.
f
V’G~°~’ = —1,
(7.6)
These equations mean that the solid angle of a closed surface is 0, —2ir or —4ir if the observation point is outside, on or inside this surface, respectively. Eqs. (3.30), (3.31) can then be written
~—fdO’ r’[v’. V’G’1’O(O’) v’~ ~—fdO’ G~’N(0’)_~_fdO’ r’[v’ V’G~’Oc(O’) —
2O~(O) 20~(O)= —
—
—
0~-
f
v’G~~’o~(o)],
(7.7)
dO’r’[v’. t7’G~’0~(O’) V’. V’G~° —
5’çb~(0)]
f
— ~—
V’~
dO’ G~’N(O’)
—
~—fdO’r’
[ii’.
V’ G~~’ 0~(0’)
—
V’~ V’ G~c’0c(O)].
(7.8)
All the integrals containing v’ V’G~B’are now regular. The two source terms containing G~~’ still exhibit a logarithmic singularity which we treat analytically. Instead of solving
f
5AB’ =
do’ G~~’N(O’),
(7.9)
we solve = 5reg + Sanai,
=
~
f
Sreg =
f
do’ ln[~2 + (0
—
do’(G~n’N(o’)+
O’)2],
~2 =
(p
ln[~
+
(0
—
01)2 ]N(o)J~
—p’)2/pp’,
where Sreg is a regular integral and Sanai is the part which can contain a logarithmic singularity.
(7.10)
(7.11)
346
/ ERA TO stability
R. Gruber et a!.
code
7.2. Discretisation In order to solve the system of coupled equations (7.7), (7.8) we expandN(0), Ø~(0)and O~(°) in terms of piecewise constant finite elements fI+i/2 Nx
N(O)
—
Nx
Nx
0~~O)= i0 ~ b
i=0 ~ a1+1/2 f1+i/2~
Ji+1/2~
1=0 ~
—
c f+l/2fJj+1/2~
(7.12)
integrate over 0 and end up with a set of coupled matrix equations 2b—2cAb—B•a—C•c,
0=D~b—Ea—F’c.
(7.13)
The matrix elements are
1
01+1
A11’= ~
1
~
2ITA0 1+112
dO
~
(f
dO fi+l/2
f
0j+1
ln[~2 +
2+
do’ fj+l/2 (G~s’+~Lln[~
0~
01+1
O’)2])
U
dO fi+l/2
(
~0i+1 dO fi+i/2
Ojf~
dO’ r v
V G~’
2ir
dO’f/+l/2 r’v’
I~0j
0
f 0
01+1
1
f
1
fl+i/2~1/
0~
61
1 2ir~01÷112 Ui
dO’ r’v’
2ir
dO’fj÷q 12r’v’ V’G~c’
=
—
I 0~j4~
dO f,+l/2
_____
CS
(0
J”
(0 _O’)2] dO1)~
6i—1
1
sc
fi÷1/28i/
—
01+1
f
B”= 2zM1~112 1 ~fi+l/2
dO’ fj÷112rv V’G~’
fi÷i/2
Oi+1
2sr
6/÷i
f
V’G~’
dO r v
fi+1/25i/ 0
—
2 + (0
do’f
—
I
I
~I~-,Q ~ ~,~S’j
~
01)2])
1÷112~G~5’+1 ln[~
~cs 0~
ln[~2+(O _0I)2] 01+1
Fe”
= 21r,~.Oj+1/2r 1
0/+1
dO fi+l/2 6
cc
do’).
I foi
27T
dO*fj+i/
f
2rIvI. V’G~~’fj+1/2~jf0 —
~.
r’v’•
V’G~~’ )~.
(7.14)
The expression v’’ V’G~B’in polar coordinates becomes
rag
aK01
~ __________
=
(a2~”~ + 4r’r)2
1 1a
2
(
2 p’cosO ,
a
ag
,
az’
op +~I
sin o
)
+ (p
2
—
p’)
+
2p[p
—
p
ap sin(O
—~
—
0)]}~
(7.15)
R. Gruber et al.
~ r, r’) = 2
a
=
=
/ ERA TO stability code
p)2 + 4pp’ sin2
—
~I
347
Note, that a can go to zero only when A and B both denote the plasma surface or both denote the wall. 7.3. Vacuum energy matrix Eliminating c from eqs. (7.13) we obtain b
a,
=
0
=
[A —21 —(C
—
21)F’D]’ [B—(C
(7.16)
21)F1 E}.
—
Whenever the wall goes close to the plasma (hEll << MAll) C~A—E,
D~A+c,
F~A,
E~B,
(7.17)
and Q becomes
(7.18)
Qn~E’B.
Note, that for a wall tight on the plasma all subdeterminants of Q are This implies that q.~0 which corresponds to the condition (3.17). The way that we calculate the matrices enables us to pick up this singular behaviour of Q and therefore it is possible with our code to put the wall as close as you like. The real and the imaginary part of the vacuum energy matrix V are defined by °°.
W~~1 = (fVX)’~”.
(7.19)
This expression must be compared with eq. (3.24)
W~’=
—~-
~
(a*i+h/2QjjaI+h/2i~0j+i/
(7.20)
2)~
and 2 @1~)1+hhI2= X,~1—
nq~-1”~i~_~I = HRX, 2 Eqs. (7.18)—(7.20) lead to VRI
=
X,
+
(aI)11h/
=
_______
—
L~x1
nq X,+ 1 +X1 2
=
H’X.
(7.21)
(7.22)
H*RIQHRI,
with —
~lJ
T
~xa+1/2
~xJ÷1/2
q2
n~ LWj+1/2
~.Wj+l/2
A~
Z4. Then—Oterms Similarly to (7.12), we expand (1/r’)v’ V’al’ into
~
1~112f
V’al’)
=~
d
1÷112,
(4w’. V’al’)
=~e1+h/2f;.+i12,
(7.24)
and obtain the matrix equations t~’Hd—Pe,
u=Rd—Se,
(7.25)
348
R. Gruber eta!.
/ ERA TO stability code
where t’ =
2
=
1
f
2sr
2ir
+
f~+i,2 27rL~0f dO
0
2ir
I
2irM
0
f f
HZ! = 2yr~0
2sr
21T
pit = ~
d01+112J —--v 0 ~dO,
f
d0’fj+112(4rr’G~5’),
f
d0’f,÷ii2(4rr’G~c’),
dO f~’+1I2 0 dOJ~+
2,r
1~2
2sr
2ir
f
~—~-~j 1 ~ dOf,+~p 0
dO’f,+112(4ir’G~1.),
0
2ir
I
v’. V’(4n~’G~5’),
2ir
2ir~O
=
~ r
2zr
f
dOf÷112
j
d0’f,+ii2(4~’G~~’).
Eliminating e in (7.25), we obtain 1R]1(r—PS~u). d [H— PS The vacuum energy (3.41) together with (7.19) becomes
(7.26) (7.27)
W~=
W~=
(7.28)
{X*(u+W+V)X]I,
with =
[~o
1~o1iJf~L~f] Sij, U
[d1dJ~o1~oJ/(~ d1~01)]o1J.
=
(7.29)
7.5. Resolution of the integrals The u integral in K~(~), eq. (7.1), can recursively be expressed in terms of elliptic integrals of the first the second (E) kind. It turns out that 4n ~+1 K0 2n+1 ~—1 2 K1=—K-i--—---—E, i~—l l—r~ —
K~+1
—
—
(K) and
2n—l K0_1, 2n+l
(7.30)
K0=K.
(7.31)
R. Gruber eta!. / ERA TO stability code
349
The n-derivative can be found by differentiating (7.24) and by using 1 ~K 1 ri + 3~ 0 1 ~E ~l E (3 + ii)Kj~ = K]. a~ 2Ø~ 1) i~ a~ 2(t~ 1) ~ —
——
—_[——
—
(7.32)
—
All the regular integrals in (7.14) are performed by means of a 4-point integration formula
f f
(7.33)
ab[f(—h~, o)+f(o, h~)÷f(h~,o)÷f(o,h~)],
f(x,y)dx dy
x—a y—b
4) which avoids any values
where = 2a~/i7~ and h~eliminating = 2b~./i7~. Eq. consists of an integration rule in 0(h on the h~ diagonal x = y, thus the (7.33) numerical singularities there. In the code, the singular part of 5afl~,eq. (7.11), is calculated analytically, i.e. 1 ~an~
=
~j+i
— —
01+1
~ a’~’~f
dOJ~÷
f
112 Ui
£
2
+
(7.34)
(0~+~O~)2]. —
dO’ 2 ln[b
Note that for the matrix B (7.14), p and p’ are on the plasma surface and at 0’ = 0, p’ = p a logarithmic singularity is obtained in (7.34). In the integral E, however, there is no singularity but for a wall near the plasma surface (when b is smaller than the mesh size) an apparent logarithmic singularity is built up numerically. For such a case it is advantageous to treat this integral analytically, to let us then solve 01+1
f Ui
o
0i+1
dO
f
do’ ln[b2
+
(0
—
01)2]
=
3(0
—
0
0~)2+ 4b(0
1+1
—
O~)arctan
Ui
2 in b2 + [—b2+ (O~+~ 0~)2]ln[b2 + (O~+~ O~)2]. (7.35) b Here b2 is considered to be constant in the interval O~~Z0 s~0~÷~ i.e. b2 = 0 for the matrix B (the value of the integral is then —(O~+~ O~)2[3 in(0 2 = (p~ Pc)2/PpPc atO = (O~+~ + 0~)/2for the matrixE. 1+1 singularities 0)2]) andinb (7.26). In the same way we treat the logarithmic —
—
+
,
—
—
—
—
7.6. Definition of the shell The position of the conducting wall is defined in subprogram SHELL. At the moment we measure the distance between the plasma surface and the shell by the quantity REXT-l, where REXT is given in the NAMELIST. The shape of the wall is defined by Pw =
Pp +
(REXT1)[p~(0)+ Pp(7T)l/2.
(7.36)
This means that we add a constant quantity in the radial direction to the plasma surface. It is evident that any other shape can be introduced by changing SHELL.,REXT > 1000 means that the conducting shell is at °° without the stabilizing line current along the toroidal axis.
8. ERATO 3 8.1. Organization ERATO 3 constructs the symmetric block structured matrices A and B which form the eigenvalue problem Ax = w2Bx. Only the upper symmetric part is stored as a one-dimensional table. One block contains all the contributions for a given s-interval. The blocks overlap due to the linear basis functions of X in the radial direction.
350
R. Gruberetal. /ERATO stability code
x,. 12 11
TT
48 47
x
46 45
x
33
44 43
x
1931
42 41
x
6~2B 40 l 15 27 39
x
2 1
l4~26 1325
x
Xl
V1V1
23 35
155
—
10 9 Xl Vi Vi X2 Y2 V2 X3 V3 ‘13 X4 V4 V4 X5 Xl Ti Bi I vi \\ vi X2 _________ B2 V2 V2 “ X3 __________ B3
________
21
8 l
I
6 5
________
1729
I
V3 X4
________ 156
1 3
0——-
“
14
B4
V4 X5
—
Fig. 6. Block structures of A and B for N~= 4.
38 37 X2
V2V2 X3
s V3V3
X4
Fig. 7. Global numbering in the
Y4V4 X5
(s, x) plane.
We obtain a matrix structure as shown in fig. 6.and the corresponding numbering of the unknowns is given in fig. 7. Two unknowns per mesh point are set for the real and the imaginary parts. In the case given we have chosen N~= 4 and N~= 5. One block has the size 8(N~+ 1) X 8(N~+ 1) which in our example corresponds to 48 X 48. Since only the upper half is stored 48 X 49/2 matrix elements define one block. The overlap goes over 2(N~+ 1) unknowns, giving 12 X 13/2 matrix elements. The length of the whole matrix which corresponds to the total number of unknowns is (3N 5 ÷1)(2N~+ 2). In the construction of the matrices we calculate all the possible contributions per mesh cell (16 X 16). Regularity conditions (only for the first block), periodicity conditions (for each block we eliminate all unknowns which lie in the domain ir <~< 2ir), and boundary conditions (for the last block only) are directly imposed on
XLj,l
~‘~‘‘/2.j.1 VL,ys ).1
4,..~ 3
2 1 ~LJ
12~16 11 15
‘1~~y2j
Fig.
8.
j’l
7
1014 9~l3 ‘~‘L.~/2,j
~i.l
0
X component
X
V and V components
0
physical components
6 ‘-‘5 XL,tj
Local numbering in one mesh cell.
R. Grubereta!. /ERATO stability code
351
these local 16 X 16 matrices. After that we add them up in a block. Two blocks are filled up consecutively. The overlapping part can then be directly added up. This is done in subroutine CONMAT. The blocks ofA are stored on disk file NDA and the blocks of B on NDB, see table 9. For the local 16 X 16 matrices containing all the contributions per mesh cell we introduce a new numbering system (see fig. 8). These numbers are related to the global ones (fig. 7) by means of the subroutine ~DEC). 8.2. Potential and kinetic energy The 16 X 16 matrices, i.e. the contribution to the potential energy matrix A per mesh cell, are built up in the subroutine (AHYBRD). In the integral (3.11)—(3 .13) we replace the vector components X, Y and V by values in a mesh cell, namely XRI
=
[X~’ + x~J÷, + X~,,~, + X~,~+11/4
axRh/ax = [X~J+
1+ X~’1,1+1 X~’ ax’~/as= [X,~ +X ,~+, X~’ —
—
—
—
2(X~+i Xi), —
X~1,1]/
x~/+i}I2(s
1+, se), —
yRI
=
VRI =
(Y~’~12,1 + Y~’~12,1÷1)/2, aYRh/ax = (Y~12j+1
—
Y~’112,,)/(x,+1
—
Xi),
(V)~12,1+ Vy~’~12,1+1)/2, av’~-’/ax=(V~’~12,1+1 V 1/2,J)/(xJ+1 —x,).
(8.1)
—
The positions of the nodal points in one mesh cell are explained in fig. 3 and the numbering of the unknowns in one mesh cell for the code in fig. 8. Note that eqs. (8.1) are identical to centred finite differences which will be applied to the integrals (3.12) of the variational form (3.11) for the mesh cell (s~~ S s1+1 ~ X~ ~
i’~ ~~-‘--+nX’~,
qax
qa~
ax’~ av” as ax I~= HXR /3XnqX’ ‘2 =
~,
jl
—
—
ax’ av’ ~,
as
ax
nq V’ ÷ i3x
~—
2 XR+a lnr2 VR ÷~
I~~a ln r as
ax
+
ax
k 2 3 4
5 7
~—,
I’~= HX’ + (3~nqXR+ nq VR + (3x
ax yR asnqY’
Table 9 IODSK3(k)
1
,
reads NEWRUN from NSAVE writes NAMELIST back on NSAVE reads vacuum contribution from NVAC writes matrix A on NDA writes matrix B on NDB reads equilibrium quantities
ax
~
as
ax
-~-
ax
+
~-,
as
(8.2)
R. Gruber et al. / ERA TO stability code
352 ,
a
a ln r2 a ln r2 ax v’+ ax
in r2
as
ax’ a v’ a “ ôs
ax
ax
i~’=x’~’. The quantity H = (q/T) a (T/q)/as
+
2j~O
2.The potential energy W,, (3.11) is then 5sr/l V~t’I
w~,= ~-ff~-dx [aI~2~J~2)
+
b(J~2~J~2~) + c(J~2+I’3~)+ d(I~2+I’)
—
e(I~2+Iç)].
(8.3)
In the same way one has to proceed for the kinetic energy, where no interactions between real and imaginary parts take place. The kinetic energy becomes
K = +ff~dx{f(XR2+X12) +g[(VR
+ yR
—
P~~)2 + (V1 + Y’ + flxx’)21 + h(YR2
+
~2)}
(8.4)
The variation is afterwards performed after all the nodal values. 8.3. Symmet7y conditions As mentioned before we only treat half the plane, and add half intervals in the symmetric half plane in order to be able to impose the symmetry conditions (3.18) which imply X’~(x)= —X’~(2ir x), yR(x) Y’~(2ir x), —
—
X’(x) = X’(2ir Y1(~)= —Y’(2ir
—
—
x)~
V’~(x)=V’~(2~r—x), V1(x)=—V’(2~i—x).
(8.5)
To fulfil (8.5) we introduce new variables XR( 2~ _X)rX’~(2iT—x)
+X’~(x)r0, x) V’~(x)= 0,
~‘(2n~ —x)
=x’(2ir —X)—X’(X) 0, X’(2z x) V’(2ir x) + V’(x) = 0,
VR(2~r x) VR(27r t(2ir—~)=Y’~(2ir—~)— Y’~(~)=0, ~“(2ir—~) Y’(2ir—~)+Y’(X)~r0, ?‘ which have to vanish. We first write eqs. (8.6) in matrix form
(8.6)
~Ux,
(8.7)
—
—
—
—
—
(x~X’~,X1, VR, V’, y’~,Y’)
and perform a matrix transformation on the eigenvalue problem Ax = w2 Bx which becomes (8.8) with
A
=
(U_i)T AlT1,
B = (U_l)TBIrS.
(8.9)
In the transformation (8.9) the block structures as well as the symmetries ofA and B are the same for A and ~&. In the code, we impose symmetry conditions cell after cell in subroutines (INTEGR) and
R. Gruber et al.
/ ERA TO stability code
353
& 4. Regularity and boundary conditions The regularity conditions (3.16) are easily imposed by X~=X’01=0,
(8.10)
0~j~N,<
throughout the subroutine (AWAY). In our code we have 2 different boundary conditions. The simple condition for a wall straight on the plasma surface (REXT < 1.0) becomes X1=X~1=0,
~
(8.11)
The other possibility is to prescribe a conducting shell at some finite distance (1.0
big value of REXT (REXT> 1000) simulates a shell at infinity. For these cases we add the vacuum energy matrix V, defined in (7.22) to the last block of the matrix A. This is done in subroutine (ADDVAC).
9. ERATO 4 9.1. Eigenvalue problem ERATO 4 solves an eigenvalue problem 2Bx (9.1) Ax’~c~., for symmetric, block structured matrices A and B. In addition, B is positive definite. The eigenvalue w2 of the problem can be positive (corresponding to an oscillatory, stable displacement), zero (the marginal situation) or Table 10 Parameters for VEKIT AL(3) AL(1) = ALO AL(2) AL(3) EPSCON EPSMAC NPIN(9) NPIN(1) = 0 NPIN(2) = N~ NPIN(3) = (3N NPIN(4) = NPIN(5) = NPIN(6) = NPIN(7) NPIN(8) = NPIN(9) =
5 + 1)(2N~+ 2) NITMAX NPRNT NDA NDB NDLT NDS
I 0 0
I I I
shifting
eigenvalue eigenvalue by normalization eigenvalue by Rayleigh quotient c.~,2 = (xT, accuracy of eigenvector machine accuracy
I
vector iteration required number of blocks
I I I I I I I
number of vector components maximum number of iteration steps printer output channel channel number for matrix A channel number for matrix B scratch unit number scratch unit and storage ofx
0 0 0 0
number of eigenvalues less than ALO number of iterations used number of converged components 0 all subdeterminants (A) regular = —1 subdeterminant (A) singular
NPOUT(4) NPOUT(1) = NPOUT(2) NPOUT(3) = NPOUT(4)
NEG NIT NCONV NSING
The eigenvectorx found at the end is stored on unit NDS and will be reread by ERATO 5.
Ax)/(xT, Bx)
354
R. Gruberetal. /ERATO stability code
negative (a purely damped or a purely growing, unstable mode). The eigenfunction x contains the displacement components arranged in the way described in fig. 7. The problem (9.1) is solved by an inverse vector iteration, namely by VEKIT, a subprogram of the block matrix library HYMNIABLOCK [10]. To have access to all the eigenvalues we perform a spectral shift. Instead of solving (9.1) we solve 2Bx, (9.2) Ax~ where A = A w~Band ~2 = ~2 ~4. By an adequate choice of w~(the NAMELIST quantity ALO) we are able to pick out each eigenvalue of the spectrum that we like. To make the iteration efficient we decompose A into —
—
A=LDLT,
(9.3)
where L is a left hand side trigonal matrix with ones on the diagonal and D is a diagonal matrix. The number of negative terms in D is identical to the number of negative eigenvalues. Therefore we always know which eigenvalue of the spectrum has been calculated.
9.2. How to call VEKIT The eigenvalue solver VEKIT is called in EIGVAL by CALL VEKIT(AL,EPSCON,EPSMAC ,NPIN ,NPOUT).
(9.4)
The meaning of the formal parameters are described in table 10 (I 10.
=
input parameter, 0
=
output parameter).
ERATO 5
10.1. Diagnostics This main program performs all the output with the exception of the eigenvalue which is printed out in VEKIT, a subprogram of the library HYMNIABLOCK [10]. We first read in the eigenfunction from disk unit NDS, separate it into the six components XR, x’, yR y’ VR and 0 which are transformed back into ~ and see eq. (3.9). These components are normalized with its highest absolute values (NORM FACTOR * 999) and the three most significant digits are printed out as a table ins and x. The s-coordinate goes from up (s = 0) to down (s = 1) and the x coordinate from the left (x = 0) to the right (x = ir). The actual values of the vector components are obtained by multiplying the tables by the corresponding NORM FACTO1~which is printed out at the top of the tables. As a next step we plug the eigenmode back into the variational form, eqs. (3.l1)_(3.15), and calculate the contribution of each mesh cell. Again we print out a table containing these cell contributions to the potential energy and to the kinetic energy. The last column at the right hand limit contains the cross sum of all x contributions for a given s interval. With such a diagnostic we are able to see where the negative driving terms of an unstable displacement come from. The test case in section 12 shows such a plot for an internal kink mode (all the negative contributions to the potential energy are internal). As an additional control, we add up the potential and the vacuum energy and form together with the total kinetic energy the Rayleigh quotient which should correspond to the eigenvalue obtained by HYMNIABLOCK. The diagnostic part concludes with a “MAP OF POTENTIAL ENERGY CONTRIBUTIONS” which enables us to find out which of the 5 integrals I~to I~contributes to the potential energy. ~,
~,
~,
~
~,
10.2. Plotting The plotter or graphics output is controlled by the quantity NFIG (in NAMELIST). NFIG = 0 no plotter output, NFIG ~ 16 displacement graphs as shown in the test run output are drawn at NFIG different angles at which the torus is cut perpendicularly to the magnetic axis (in the s—x plane),
R. Gruber et al.
/ ERA TO stability code
355
ARROW enables the size of the arrows to be changed, ANGLE (NFIG) denotes the angles in degrees at which the NFIG torus cuts have to be drawn. ANGLE is given in NAMELIST. The quantity ITEST is CDC dependent. It must be set to 1. ITEST = 0 only checks if the plotter output is in the limits without drawing anything. The plotter or graphics output passes through the following subroutines ABSDIR calculates the directions and the amplitudes of the displacement in each mesh cell, PLOTER organizes the plot, FLECHE draws an arrow, SURF draws the plasma surface, TITRAX produces titles and axes. The plotter or graphics output contains three machine dependent subroutines which must be adapted for each computer: PLOTS initializes the drawings, PLOT moves the pen up or down, SYMBOL draws a string of characters. Details about these three subroutines are given in the listing of ERATO 5. To avoid using another machine dependent subroutine (NUMBER) digits are converted into characters by ENCODE which is not standard FORTRAN but available on many computers. The character obtained is then drawn by SYMBOL. The CDC subroutine DATE is also used.
11. Test cases 11.1. Preset values (internal mode) The code can be run without giving any data between the $NEWRUN and the $END card. This default case calculates an internal mode of the Solovev equilibrium, eq. (2.8), with the parameters given in table 5 of section 6. As a result one obtains an eigenvalue w2 = —1.5883. Note that the $NEWRUN card must be preceded by 4 label cards. Setting NTCASE = 1 in the namelist gives the result w2 = —1.5889.
11.2. Test run output (kink mode) As a second test case we treat a free boundary (REXT = °°)kink mode with a growth rate _~2 = 1.2373. This case is shown as test run output at the end of the paper. Only the first two and the last two pages are shown. As input card stream the card sequence 1st card: this is the ERATO test case 2, 2nd card: an external kink mode (n = 2, E 2 is calculated), 3rd card: note that these 4 cards must always 4th card: precede the NEWRUN input, 5th card: $NEWRUN, 6th card: NPSI = 14, NCHI = 15, WNTORE = —2., 7th card: ELLIPT 4., NTCASE 2, ALO —1.2, 8th card: QIAXE 1.6667,REXT= 1001.0,NMESH= 1, 9th card: NFIG = 4, ANGLE(1) 0., 30., 60., 90., 10th card: $END must be given to obtain the test run output. For this case we also show the plotter output of the resulting typical kink-like eigenvector (fig. 9). The (~,~) displacements are shown at each (s, x) mesh cell centre at (NFIG = 4) differentq angles (ANGLE(l)= 00, 300, 60°,90°).Changing REXT = 2.0 gives an eigenvalue w2 = 1.1015. I~2
356
R. Gruber eta!.
~ /1 ~ IT~V.’,f~~
code
/~\‘-/~\
ii
/
--
/ ERA TO stability
‘
~ IV
/‘~,,/:
~
/~st~
~
~
~
~
:
~‘
~ço ~
:::::: --
~
~:-
~ I
~
LPUSPNNE
-
/L
~
~
V
if
/
- -
~~
“,
I
-
L6US6NNE
.
S
‘
,
—
.•‘~
:~O~j , c~—~:~_
5
\‘
— —
—.~
~
\~ ~H~~•’’’-\L’ ~
1
\V~~r
/
-
/
/
/
f
~EP6TO
=
/
/
~j~=-
60
26/03/86
/
-
/
\V
P1-/I
—
-
_\
\
~
rpr
26/03/80
—\
——S
// / / 1y~~//~ ~ / V \l~ .-‘~1~’/ ‘I ,‘~L~1~1 ‘I \~ S~,y
~l4
30
-N-\
~
\
/1~
I
‘ERPTO
~
~~‘-~‘‘
I
LRUSRN~IE
—
~ /~ ~~ I’H~~\
V I~~~V_: ~
CPPP
-~
IV
_—_/7~
~
\V
EPF
CPPP
~\
~
-‘ -
-
26/03/80
/
-
~‘-.-,,
1
V
/ -
:~~
PHI
‘S
L’-
H
ERPTU~
~
\
=
‘
7.I~
\
‘~ ft
I V~’~ ~
£PF
,/
~“ -
‘~
_
p~
CRPP
~
-
:~ ‘
\~
~
CRPP
£PF
—
LRUSRNNE
PHt
‘cR610
Fig. 9. Plotter output showing the displacement vector field at different ~ angles.
=
90
26/03/86
R. Gruber et al.
/ ERA TO stability code
357
All these cases are related to the calculations performed in ref. [27]. The computing time of such a case is around 2 mm on our Cyber 170-720.
12. Instructions to the user
12.1. Convergence properties A regular finite element approach [20] always shows quadratic convergence behaviour when the number of intervals goes to infinity. As presented in a previous publication [8] we have to take a very large number of mesh cells to obtain the required quadratic convergence and, sometimes, the actual computers are not big enough to solve the stability problem by such a method. This inconvenience comes from the fact that the eigenvalue is always approached from the stable side of the spectrum, i.e. from above. A coupling between the unstable, not yet enough converged mode, with stable, almost marginal modes of the continuous spectrum, can strongly affect the convergence behaviour. In the finite hybrid element method, however, we reach a lower energy level, and therefore prevent the coupling with the continuum. But another problem arises which we call the spectral shift in q [15]. However, the elimination of this spectral shift, eq. (4.17), is only valid in the case of a Solovev equilibrium. We have seen that this shift depends on the shear at the plasma surface. A more general expression has still to be found. In fixed boundary cases we even sometimes find higher order convergence laws. Such a higher order convergence is only found in the absence of the vacuum region since the vacuum contribution in its actual form always shows a quadratic convergence law. For Tokamaks it is seen that the maximum F2’s are always higher in the discrete case than for the exact value. The marginal situation of the continuous spectrum can also become slightly unstable and a convergence study must push this apparent unstable mode back to the stable region. A careful study of the localization of such a mode by means of the potential energy and the kinetic energy tables can be a helpful tool to recognise a continuous mode in the unstable region. As a consequence, the user must be very careful in the interpretation of his results. 12.2. How to switch from equilibrium to ERA TO As mentioned before, ERATO needs to know t~(r,z), p(i,li), T(’I’) and TT’(i~),such that the equilibrium equation = —r2p’ TT’ (12.1) —
is fulfilled. Whenever the equilibrium code is not written in this “natural” system [18], one has to perform a transformation of p’ and TT’ in the subprogram TANDP. The normalization with respect to the magnetic axis R and the normalization of T with respect to the magnetic flux at the axis T 0 is done in subprogram NOREPT -and the user has nothing to care about. As a default value, we use unit NDORY = 11 (catalogued on local file name DORY) to pass the necessary quantities from the equilibrium code to ERATO. The accuracy of a numerical equilibrium using a finite difference method is somehow related to the accuracy of a finite element method as used in ERVATO. As a rule of thumb this relation is essentially given by Nr
6N5,
N2
3M5.
(12.2)
These relations come from the fact that in a finite difference approach the derivatives are found in between two intervals, whereas in the finite element method only one interval is used. In addition the number of r-intervals must be doubled in respect to the number of z-intervals, since only the upper part of the torus is treated. Another
358
R. Gruber et a!.
/ ERA TO stabi!ity
code
factor of 3/2 is introduced to take into account that on the average about 1/3 of all the r—z mesh cells are outside the plasma region. Note that in the high i3 case where the magnetic axis is pushed outwards the factors in (12.2) must be increased. 12.3. Choice of NAMEL1ST parameters All the NAMELIST parameters can be modified in all the 5 main programs. Generally one defines them in ERATO 1 and leaves them unchanged throughout the whole tr~atment.However, it can sometimes be of some interest to build the matrices A and B only once and go through the spectrum with different values of ALO c4. The matrices A and B are then stored on disk files NDA = 7 (with local file name TA) and NDB = 9 (TB) and the NAMELIST quantities are catalogued on unit NSAVE = 8 with the local file name OUT. In the NAMELIST input, there are some special points to be considered: (I) Convergence studies are simpler when the same number of intervals in s and x is chosen, i.e. when N~=
N5
+ 1. (II) For an unstable configuration ALO <0. (III) A fixed boundary case yields REXT < 1; REXT> 1000 means a wall at infinity and in between there is a wall at a distance (REXT 1) X (plasma diameter at z = 0)12. (IV) The quantity EPSMAC has to be set to the machine accuracy. For a CDC 6500, EPSMAC = 10_12. (V) The mass density p is preset to i in the whole plasma region. When a mass density profile p(s) is required one has to foresee this in the subprograms EQANAL and NUMEQU when filling up EQ(7,JCHI). The user —
EPF LAUSANNE MOS/BE 1.3 1EV 48~—0621f07180 12.12.28.72561gM FROM 12.12.28.IP 00071232 WORDS — FILE INPUT DO QL 12.12.28.T2551,CMILi0000,Ti000. GRUBER,I~RPP 12.12.28. 1 12.I2.30..REQUEST,SOUROE,~PF. 12.12.30.UPOATE.F.N=OLDPL,C=0,L=0,S. P i2.12.30.REPARES OLDPL 2 12.12. 31..UPDATE OREATION RUN 12.i2.3i.CREATIN~ NEW PROGRAM LIBRARY 12.12.Li9. UPDATE COMPLETE. 12.12.49.CATALOG, SOURCE,ERATO, IO=LISRS. PP=999. 12.12.’.3.INITIAL CATALOG 12.12.51.CT 10= LIARS PFN=ERATO 12.12.51.CT CY= 001 00020 PBS 12.12.5I.UPOATE,Q,L=Q. R 12.12.51.OUTINES FOP ERATOI 3 l2.~3.10. UPDATE COMPLETE. I2.~I3.I0.FTN,I=COMPILE,L=0. 12.13.10. Li 12.i’..kS. 32.587 OP SECONDS COMPILATION TIME 12.14.49.MAP(PART) 12.1~.Li5. 5 12.14.le.LGO. ~ ERATO I ~ 12.1Li.1.5. S 12.15.22. STOP 12.15.22.0554000 CHAMP MEMOIRE MAXIMUM 12.15.22. 27.007 OP SECONDS EXECUTION TIME 12.j5.22.RETURN,LGO. 12.15.22. 7 12.15.22.UPOATE.9,L0, R 12.i5.22.OUTINES FOR ERATO2 8 12.15.28. UPDATE COMPLETE. 12.15.28.FTN,ICOMPILE,L0. 12.15.28. 9 12.16.08. 12.095 OP SECONDS COMPILATION TIME 12.I6.08.LGO. ~~4” ERATD 2 12.16.08. 10 12.16.13. STOP 12.16.13.03560/lB CHAMP MEMOIRE MAXIMUM 12.16.13. .145 OP SECONDS EXECUTION TIME 12.16.1 3.REWIND,LGO. 12.16.13. ii i2.i6.13.UPOATE.9,L0. P 12.15.13.OUTINES FOR ERATOS 12 12.16.19. UPDATE COMPLETE. 12.16.19.FTN,I=COMPILE,L=0. 12.16.19. 13 12.16.54. 11.887 CP SECONDS COMPILATION TIME 12.16.54.LGO. ERATO 3 12.16.54. jL 12.17.22. STOP Fig.
12.17.22.OT6k008 CHAMP MEMOIRE MAXIMUM 12.17.22. 10.880 CP SECONDS EXECUTION TIME 12.17.22.REWINO,LGO. 12.17.22. 15 12.t7.22.UPOATE,Q,L0. P 12.i7.22.OIJTINES FOR ERATO4 16 12.17.29. uPDATE COMPLETE. j2.17.29.FTN,I=COMPILE,1=0. 12.17.29. 17 12.18.06. 6.362 OP SECONDS COMPILATION TIME 12.18.06.FILE.TAPEIO.RTS. 12. 13.06. 18 12.18.06.LDSET.FILESTAPEIO. 12.18.06. 19 12.18.06.LGO. ~ ERATO L 12.18.06. 20 12.22.36. STOP 12.22.36. 86.5L+L. CP SECONDS EXECUTIOP~TIME 12.22.36.i0Lii000 CHAMP MEMOIRE MAXIMUM i2.22.36.REWIND.LGO. 12.22.36. 21 12.22.36.UPDATE,9 1=0. P 12.22.36.OUTINES ~0R ERATO5 22 12.22.Li7. UPDATE COMPLETE. 12.22.Li7.FTN. I=COMPILE,L=0. 12.22.47. 23 12.23.19. 14.252 OP SECONDS COMPILATION TIME tZ.23.l9.LOSET,FILESTAPEIO. 12.23.19. 24 12.23.19.IGO. (RAT/I 5 12.23.19. 25 12.23.24. NOil~FATALLOADER ERRORS — SEE MAP 12.23.32. STOP 12.23.32.0644008 CHAMP MEMOIRE MAXIMUM 12.23.32. 2.604 OP SECONDS EXECUTION TIME 12.23.32.EXIT. 12.23.32. 26 12.23.32.QP 00010816 WORDS — FILE OU1PUT, DC Lie 12.23.32.MS 86016 WORDS ( 526848 MAX USED) 12.23.32.CPA i84.q58 SEC. 6.165 AD.). 12.23.32.CPB 59.470 SEC. 1.982 ADJ. 12.23.32.10 160.905 SEC. 10.726 80.). 12.23.32.C$ 1075i..317 045. 388.972 AOJ. 12.23.32.SS~-~COIJT OU CALCUL 285.90 FR ~‘ 12.23.32.00 PERIPHERTE EN SUS (MESS CC. TRACEUR) 12.23.32.PP 308.976 SEC. DATE 08/09~80 12.23.32.CC 25.70 FR I/O 0 PAUSES 12.23.32.EJ END OF JOB. ~‘ rn~8~4*~#
10. Control stream of ERATO.
T2SSI9M 125519M
I//f If/f
END OF LIST I//I END OF LIST If/I
R. Gruber eta!.
(VI)
(VII) (VIII)
(IX) (X)
/ ERA TO stability code
359
is warned against too realistic mass density profiles which fall to zero at the plasma boundary. These can cause unrealistic values of w2 which tend to infinity with 1/p (see ref. [19]). NLY = T calculates the whole problem, NLY = F eliminates the toroidal component Y by imposing V• = 0 in the potential energy and Y = 0 in in the kinetic energy, NLTORE = T for toroidal geometry (tokamaks, spheromaks, reversed field pinches), NLTORE = F for helical geolnetry. This case is not yet included in the published version of the code. NMESH = 0 for an equidistant s-mesh, NMESH = I for a non-equidistant s-mesh. Accumulation around integer nq surfaces is made. The user can easily change the accumulation algorithm in MESH. WIDTH 00 for equidistant x-mesh, WIDTH = finite accumulates the x-mesh at the outer region. The user can again easily change this. 1—0 operations in IODSK3 of ERATO 3 and in IODSK4 of ERATO 4 should be accelerated by introducing direct access of the matrix block records and by directly buffering into main memory.
10.01
ERATO1
1.1.02
CLEAR
(cocJ LOCE
11.03
PRESET 1.1.01
LABRUN 1.1.04
DATA 1.1.05
AUXVAL 1.2,01
ORGNIZ 2.0.01
ERATO2 3.0.01
ERATO3 4.001
ERATO4 5.0.01
ERATO5
2.2.01
VACUUM 3.2 01
AANDB 4 2.01
EIGVAL 5.2.01
SORTIE
Fig. 11. Overall structure of ERATO. The main programs ERATO2—ERATO5 are dummies, and the subroutines VACUUM— SORTIE can all be included in ERATO1 if sufficient main store is available.
360
R. Gruber eta!.
/ ERA TO stability code
Table 11 Adjustable COMMON dimensions No.
Block name
Variables
Symbolic dimensions
2.6
COMVEC
XT, UT, VT
(6N
X,U,VA,VB,C A1-A6 BETRZ, CHIRZ RRO, RRS, RZS, TETA
NRZ
131
lNr,Nz] Nr
(65,33)
9.1
COMMAP
CPSEQ
EQPS,EQR,EQT,EQZ ILOC 9.2
COMVID
8(Nx+l) 5 + 2)(Nx
+
1)
Default a) 1376 128
65
15
Nx
G A—F WVAC H
[2Ivx
[4Nx_2,4Nx_2] (2Nx 2)(2Nx 2)
— 2, 2Nx — (2~’~T~ — 2, 2Nx —
—
2,61 2]
—
9.3
COMMTR
A, B
4*(Nx + l)*(8Nx + 9) + 1
9.5
COMPLO
CNRRZR X XS,YS KK~KW~RK~RW}
ENs+1,Nx+l]
DW
[Ns,Nx+l,61
SUM
[N
(6Ns+2)*(Nx+l) 2*NRZ ENs+1,Nx+1J
5+1,6}
(28, 28, 6) (28, 28) (58,58)
784 8257 (15,16)
1376 262 (15,16) (14,16,6) (15,6)
a)Dej~auItforN~.:65,N~=33N~ l4,N~~’ 15,NRZ= 131.
12.4. Control stream and data set As an example (fig. 10) we show the control stream of the Solovev test case. All the 5 ERATO main programs are corrected by a batch editor UPDATE, compiled and executed. The control cards preceding ERATO 4 and ERATO 5, i.e. the FILE and the LDSET cards are set to eliminate problems with the BACKSPACE command in IODSK4. After the control stream, the UPDATE modifications for all the 5 main programs follow between the first and the second correction decks the NAMELIST preceded by the 4 label cards is inserted. They are read in ERATO 1. All these card decks are separated by end-of-record cards.
12.5. Dimensions The dimensions which have to be arranged are given in table 11. In the original version the default values are = 14,N~= 15,J\1r65 andN5 =33. These arrays have also to be modified when bigger cases should be treated. For the surface fit we have foreseen NRZ ~60. N5
R. Gruber et aL
/ ERA TO stability code
361
Table 12 ERATO 1, subprogram names and identification Name
No.
Number of
Contents
arguments
LABRUN CLEAR PRESET DATA A(JXVAL ORGNIZ FORMEQ MESH KERTP KERTES SCALES FITAXE NOREPT TANDP SELECT ISTPSI EQANAL ALPHA
1.1.01 1.1.02 1.1.03 1.1.04 1.1.05 1.2.01 1.2.02 1.2.03 1.2.04 1.2.05 1.2.06 1.2.07 1.2.08 1.2.09 1.2.10 1.2.11 1.2.12 1.2.13
PSRHO
1.2.14
5
labels the run clears variables and arrays set default values defines data specific to run set auxiliary values organizes the calculation organizes equilibrium mapping set up mesh prepares for NTCASE = 1 (Solovev analytic) prepares for NTCASE = 2 (Solovev numerical) scales the numerical equilibrium fits equilibrium around axis normalization to ERATO standard set T(~~), p(9) and p(p) selects analytic or numerical equilibrium radial counter of (v,, z = 0) analytic equilibrium quantities calculates q and non-orthogonality givesrandz(s, 8) ~p(r,z) by fit r, z derivatives of 1~/(r,z) interfaces numerical equilibrium — ERATO advances on s = constant line interpolates on s = constant
1
set up numerical EQ
1
PSIRZ
1.2.15
1 2 1 2 5 2
DERIVE
1.2.16
2
INTFAC
1.2.17
1
ADVANC INTPOL NUMEQU
1.2.18 1.2.19 1.2.20
5
SURFIT TRIBAC MERCIR INTGRA QUADRA
1.2.2 1 1.2.23 1.2.24 1.2.25 1.2.26
BALOON ALDLT
1.2.27 1.2.28
5
BETAST
1.2.29
1
DOMAIN SMLEQU
1.2.30 1.2.31
YYP BNDRY
1.2.32 1.2.33
7
OUTPUT IODSK1
1.3.01 1.3.02
1 1
6 6 3
s, x
fits the surface solves Ax = b Mercier’s criterion calculates integrals for Mercier criterion polynomial fit for MERCIR
ballooning criterion decomposes A = LDLT beta values domain for analytic fit
analytic equilibrium in domain parabolic fit dummy routine prints output disk operations for ERATO 1
12.6. Space and time requirements The memory space of ERATO can be estimated by the formula
W~l6OOO+(N~+l)(32iVx+ 18N8+42).
(12.3)
This yields 49000 60-bit words on a CDC 6500 for a 24 X 25 case. The computing time is given by T~64N~Ns(tiN~ + t2N1~+
t3),
(12.4)
R. Gruber eta!. / ERA TO stability code
362
11213
~LJ~ J1224
2O3
Jl
I
1225
1~.1~C!R
MESH
INTGRA
Jl.2.04
I
Ji
1210
I
SELECT
11302
]
7
J1.2 05 I KERTES
302
1201
ORGNIZ
1252 FORMEQ
1219
IODSKI
INTPOL
]
1 3.02
BND~J
IODSK1
1D~IVE
1
ISTPSI
J1218
7
I
Ef7
1
PSRHO 11216
I
11215
fl~7PSIRZ
7 DERIVE
1 2 29
I 7
~‘l
ADVAI’C
[1214
BETAST
11302
I 7
IODSIO
II 223 7~TR1BAC I 11.219
7
INTPOL
1211
DOMAIN
ISTPSI
1207
7
IODSK1
H
VYP
302
1,233
J1230
__________~
232
NUMEQU
ALPHA 1126
QUADRA
7
220
7
~
[i~211
[j~~J~
KERTP
7
1226
~1
I
1u~~
T 1
M3REPT
l12~
12.15
7 SMLEQU
_________
J1206 ~
I
SCALES
_Ji
2.06
I_________ 112
21
~‘[
LSURFIT
J7
]12.17 ~ INTFAC
]13 01
27
L
112 28 ALDLT
1223 TRIBAC 1219
BETAST
_________
Fig.
where ~
ISTPSI
ADVANC
11229
7
211
1 2.18
7 OUTPUT
7 SALOON
Ji
12. Subroutine map of ERATO1.
is the number of iterations in the eigenvalue solver. The time constants are t
r2
1
23 .ts when machine
coding is used for the matrix decomposition, 50 jts and t3 150 is on a CDC 6500. For a iest case with 24 X 25 intervals, the decomposition takes 9 mm, the 4 iteration steps take 3~min. All the rest (i.e. mapping,
vacuum, matrix construction, shift and diagnostics) take 2~-min. Altogether, such a case takes 15 min on a CDC 6500 which corresponds to about 1 mm cpu time on a CDC 7600.
R. Gruber eta!. / ERA TO stability code 2.3.01
363
MOOS
IODSK2 2.2.02
CALD j22O4
TETMSH 2.2.03
2.2.06
ABCDEF
[__ji~
ROTETA {
CONCEL
SHELL
] r~i
2.2.12
EKN
COC
[2.2.0) VACUUM
LOCF 2.2 07
2.2.08
QCON 2.2.09
2.2.14
WCON 2.210
IMOC 2.2.11
SYMETR Fig.
ZERQJ
MULT
13. Subroutine map of ERATO2.
Table 13 ERATO 2, subprogram names and identification
Name
No.
Number of
Contents
arguments VACUUM TETMSH ABCDEF ROTETA SHELL CONCEL QCON ZERO
2.2.01 2.2.02 2.2.03 2.2.04 2.2.05 2.2.06 2.2.07 2.2.08
WCON
2.2.09
SYMETR
2.2.10
MULT
2.2.11
4
product of two matrices
EKN EK IMGC LODSK2
2.2.12 2.2.13 2.2.14 2.3.0 1
4 3 4 1
elliptic integrals completes elliptic integrals inverts a matrix
organizes vacuum calculation
3 5 3 1
8-mesh in vacuum set up matrices it, B, C, D, E, F prepares p. 0 fdl~integration defines position of the shell contribution per cell
calculates Q-matrix reduces for n = 0 total vacuum energy matrix imposes symmetry conditions
disk operations for ERATO 2
[
1
EK
R. Gruber eta!. / ERA TO
364
stability code 3.2.10
DEC 3.2.04
AHYBRD 3.3.01
3.3.01
1005K3
3.2.01
3.2.02
AANDB
100SK3
3.2.03
CCt’JMAT
INTEGR
3.2.05
BHYBRD 3.2.07
TRANSF 3.2.09
AWAY 3.2.06
STORE 3.2.08
A~VAC Fig.
33.01
IODSK3
14. Subroutine map of ERATO3.
Table 14 ERATO 3, subprogram names and identification
Name
No.
AANDB CONMAT INTEGR AHYBRD BHYBRD STORE
3.2.01 3.2.02 3.2.03 3.2.04 3.2.05 3.2.06
TRANSF
3.2.07
ADDVAC DEC IODSK3
3.2.09 3.2.10 3.3.0 1
Number of arguments
3 2 2 2 5 3 3 1
Contents
organizes matrices A and B creates blocks of A andB cell contribution and conditions 16 X 16 function forA 16 X 16 function for B adds cell contribution transformation ofA and B removes a column and a row address correspondence disk operations for ERATO 3
Table 15 ERATO 4, subprogram naitses and identification Name
No.
EIGVAL SET3
4.2.01 4.2.02
Number of arguments
Contents
organizes eigenvalue calculation set quantities for VEKIT
3.2.09
AWAY
R. Gruber eta!. / ERA TO stability code
365
MOO?
FXFU MDO3
M012
IODSK4
BSI-UFT MDO4
MOOS
CONALR
CALD
MD12
MD12
IODSK4 ~
4.2.02
MDO6
5ET3 j4.2.O1
IODSK4 MOb
CBXLU
GETRG
MOO)
EIGVAL
M011
VDUT
PUTRO M007
CDLHXV
M012 .
IODSK4 MO1O
GETRG M008
MOll
CONORM M009
PUTRG M012
EIGEN
IODSK4 MOlO
GETRG Fig. 15. Subroutine map of ERATO4. Table 16 ERATO 5, subprogram names and identification Name
No.
SORTIE DEBUT BACKTR ABSDIR VACPOT VALUES DIAGNO PLOTER
5.2.01 5.2.02 5.2.03 5.2.04 5.2.06 5.3.01 5.3.02 5.3.03
FLECHE SURF
5.3.04 5.3.05
TITRAX
5.3.06
Number of arguments
Contents
organizes printing and plotting
prepares vector components back to real displacements prepares arrows
1
3
calculates vacuum contribution prints value table diagnoses energies organizes plot plots an arrow plots the plasma surface plots title and axes
366
R.
Gruber eta!. / ERA TO stability code 5.2.02
3.2.10
DEBUT 5.2.03
BACKTR 5.3.01
VALUES 5.2.01
5.2.04
SORTIE
ABSDIR 5.3.03
PLOTER
DEC 7
I
(GRAPHICS] .
PLOTS [GRAPHICS]
PLOT 5.3 04
FLECHE 5.3.05
SURF 5.3.06
TITRAX 5.3.02
DIAGNO
5 2.06
VACPOT
~RAPHICS]
PLOT (GRAPHICS]
PLOT (CDC]
DATE [GRAPHICS]
SYMBOL
Fig. 16. Subroutine map of ERATO5.
Acknowledgements We would like to deeply acknowledge the continuous support by the applied mathematics group at EPF Lausanne, especially Prof. J. Descloux, Drs. H. Evéquoz, Y. Jaccard and J. Rappaz. Fruitful discussions with Drs. T. Takeda and T. Tsunematsu brought one author (R.G.) to the correct formulation of the finite hybrid elements. A strong competition with the Princeton MHD group forced us to do our best. Many helpful propositions and critics were given by Drs. K. Appert,~L.A.Charlton, R.A. Dory and J.L. Soulé. All this work would not have been possible without the excellent service given by the Computing Center at EPF Lausanne. This work was supported by the Ecole Polytechnique Fédérale de Lausanne, The Swiss National Science Foundation and by Euratom.
References [1] A. Sykes and J. Wesson, Nucl. Fusion 14 (1974) 645. [2] G. Bateman, W. Schneider and W. Grossman, Nuci. Fusion 14 (1974) 669. [3] M.S. Chance, J.M. Greene, R.C. Grimm, J.L. Johnson, J. Manickam, W. Kerner, D. Berger, L.C. Bernard, R. Gruber and F. Troyon, J. Comput. Phys. 28 (1978) 1. [4] I.B. Bernstein, E.A. Frieman, M.O. Kruskal and R.M. Kuisrud, Proc. Roy. Soc. (London) Ser. A244 (1968). [5] K. Appert, D. Berger, R. Gruber and J. Rappaz, J. Comput. Phys. 18 (1975) 284. J. Rappaz, Numer. Math. 28 (1975) 15. [6] R.C. Grimm, J.M. Greene and J.L. Johnson, in: Method of computational physics, vol. 16 (Academic Press, New York, 1975) ch. 4. [71 R. Gruber, 1. Comput. Phys. 26 (1978) 379. [8] 0. Berger, R. Gruher and F. Troyon, Comput. Phys. Commun. 11(1976) 313. [9] L.S. Solovev, Zh. Eksp. Teor. Fix. 53 (1967) 626; JETP 26 (1968) 400.
R. Gruber et al. / ERA TO stability code
367
[10] R. Gruber, Comput. Phys. Commun. 20 (1980) 421. [111 J.D. Callen and R.A. Dory, Phys. Fluids 15 (1972) 1523. [12] B.B. Kadomtsev, Reviews of Plasma Physics, Vol. 2 (Consultants Bureau, New York, 1966) p. 155. 1131 E. Rebhan and A. Salat, Nuci. Fusion 16 (1976) 805. [14] K. Appert, D. Berger, R. Gruber, F. Troyon and K.V. Roberts, Comput. Phys. Commun. 10 (1975) 11. [15] R. Gruber and F. Troyon, Proc. IRIA, Paris, eds. R. Glowinski and I.L. Lions (Springer, Heidelberg, 1979); Lecture Notes Phys. 91(1979) 288. [16] M.H. Hughes, K.V. Roberts and G.G. Lister, Comput. Phys. Commun. 10 (1975) 167. [17] D. Dobrott, D.B. Nelson, J.M. Greene, A.H. Glasser, M.S. Chance and E.A. Frieman, Phys. Rev. Lett. 39 (1977) 943. J.W. Connor, R.J. Hastie and J.B. Taylor, Phys. Rev. Lett. 40 (1978) 396. [18] E.S. Weibel, Am. J. Phys. 36 (12) (1968). [19] K. Appert, D. Berger and R. Gruber, Phys. Lett. 46A (1974) 339. [20] O.C. Zienkiewicz, The finite element method, 3rd ed. (McGraw Hill, London, 1977). [21] S. Rousset, R. Gruber and F. Troyon, ZAMP 30(1979) 862. [221 F. Troyon, L.C. Bernard and R. Gruber,Comput. Phys. Commun. 19 (1980) 161. [23] R.G. Bateman and Y.-K.M. Peng, Phys. Rev. Lett. 38 (1977) 829. [24] K. Lackner, private communication. [25] T. Takeda and T. Tsunematsu, to be published. [26] C. Mercier and H. Luc, The MHD approach to the problem ofplasma confinement in closed magnetic configurations, Lectures in Plasma Physics (1974), EUR 5127e. [27] D. Berger, L.C. Bernard, R. Gruber and F. Troyon, ZAMP 31(1980)113.
368
R. Gruber eta!.
/ ERA TO stability
code
TEST RUN OUTPUT THIS IS rr1+ AN
k 4 1
3
I9~,+~2I
IXTISNAL KIN~( Iluti
NOTE
TruSt
THAT
PRECtuE Tic —.17654*06
H
TIlT CASo IS LALCILATIU
4
LAI2S
lUST AL,AYS
A
S
1 S I
I
iNPuT
_.6U9A1E~Li —.104791*0J
—..,?01r,t—ruI .I30354*i0
—.3545c+—0Z —.74s33u.—.~l. .112L’44
—.LluO{)E*40
—.c132?t*u4
—.17504,..Uu
—.184881*LrU —.I54941*uO —.69tj~4E—rJI
—.Id4A4L+.r~ —. 9Ji*uu —.bi*,i—uI
—.10113.040 —. I ISocIc —.3+4101—41
4
26
.3,7b41.~r.r2
31
—.10940+—ui —.1ru1,6L—UI —.213531—0+—. 07905+—ui .L bb+’u., .141041+lu —.I7043~~u* —.isu~,.cuu —. *.riL2o...JJ —.+7*91.uu —.i391l 1ou—.~. 27u,I+oj —. 23531+—ni —.013751—nil
—.254t~.L—r.r~. —.944601—ui —.lJrr*jri —.i021c*uu —.I7042+*u0 —.II544c~uu .o1415c—U2
—.326i71—uL —.1or3~3E*0u —.IS~r.u01*uJU —.1d335**UU —.473~3L*44 —.Iu9,+5+*Uu
—.397731—UI —.lU7Ibl*OU —.159221*00 —.i84301+00
—.15950**r30 —.S9LSbE—01
—.468811—lA —.11313c*UU —.163161,00 —.185U51*urJ —.165221*06 —.884611—Ui
~.53930lr01 —.II9U5E~uU —.1.66831*04 —.L83i.8E~U0 —.16037E*U*r —.770191—01
is
40
FIr
I
2
nn.
4. i84Eu1
9.8241—11
in.
no.562+—uI
0.5101—ui. 7.6431—UI 3.7754—UI
b.b24E—OA 7.7561—01. 8.8891—61
9.793+—Li
3.9081—UI
L.uOZL*UU
i.i33c+nnnn 1. 700+ *00
+.1041*U0 £ • 2171 *~JU
1.4151*00 1.2291 *40
5.7371—UI 7. 869s—UI 9.0021n—U1 1.0131*00 1.1271*60 1. 2401 *00
4.1791—LU
9.8241—Il
U.
1.d52i—Lni U.
4. 4.16?+—us
3.3 j.
3. 4.1311—10
9.l07cLn3. 5.
5. 718+—uI 6.6501—li l.4631—sI 9.115+—nil .nnni 1. i3njc *uU 1.2Slni*oo
5. 5,~+—~i o.55.i+-o.,1 i.sOOo—si 3.4.?OL—C1 4.13+101,0 I. 149*. •so 1.4o3,.*,.,)
rn.9441—sl. 7.u? 14—ui 0.213+—ok 9. 3’,2.—niA A.o47n.*jo I • iou •n.’I i.2741*ni4
6.457+—ui 7. list—uI 3.3215—ui H.450*.—sI 1.0 *44
o.illt—u+
6.194+—IL
5.
1. jsJs—UL
?.4L5+—uA. 5.549+—nI 9.5311 —ni 4..n..30 1 • 133+ 045
7. 530*—UI
-no.333+—ui.
3. 4.1911—Is
—i.442*n—1
4.
Cm0
IAXIS
1. A7n~.030 i.~h)t*~u
2.4~0c—ni 3.Subl—Ji 4.574+030 A • £3 3,. *‘~o +.cO?L+oi
397n.—oi
39
F IT 4.it.7+—ui U.
u. 4.167+—UI
0.
4.407+—u 0.
n.~ 0.
...
—1.442+—,.4
FIT o. 0.
I. 4541—11
4.1041—UI
I.409E—i0
~.
o.27in.—nnni 7.iulc—ni1 3.4*6+—Il 9.1651—ui 1.4595*10 i.ulu+*uni i.1741*uAi.inOj++lu u.2~5i*uu l.2,7++ou
6.2440—41 ?.4101—uI 5.5494—UI 9.+0i*—~i i.441104U 1.1931*04
6. 3971—niL 7.530+—ui 8.6+21—Ui 9.715101 1.4931+00 1.2561+44
no.5101—ui 7.6431—Ui no.775+—UI 9.9081—uI i.1U41*U0 1.2171*44
6.6241—01
6.7371—UI
7.7561—UI 4.8811—Ui 1.UULLUUU 1.1154*40 1.2291*uo
7.4691—UI 9.0021—li 1.0131*60 1.1271*40 1.2401*00
1. 3.3,1+—is
0. 0.
1.454+—il
3.4?Si—iU
1.4091—14
0.
ni.
.,.,441—u,
4.1875—si
3.
U. 3.3’ist—Io
5.710+—Ui 5.450+—on. 7.963i—ni1 9.1151—UI 1.025u*44 1.i38c*uu 1.1511*00 39
5.031+—nil
5.,4+—s1 7.377+—SI o.2491—ol 1. 342+—si
o.u37s—Ui 7.LSr.nL—nnl 4. 32+1—ui 4.4o5i—ui
4.46l+—Ui
u. 4.ls7n.—~I
—1.o521—u1 U.
‘.
is 71—ni
0.
148
14015.
o.9#nj+—ll u.0561—i1 5.Ziot.—ui i.435+’sU I.i49co~J 1.253+003
I.347+.s0 1.1614*54 i.274~*1u
FIT U.
—no.3332—uI 0.
4.15714+ 0.
R, GruberetaL /ERATO NIFAC NTFAC
NIFAC OTFAC rITFAC
OTFAC 9TFAC nIFAC *ITFAI.
NTFAL OTFAC NIFAC ,*TFAC
NTFAC rIrFAC NTFAC
rITFAC N1F4C NTFAC $5181
SITFAC NTFAC
uTP&C
Nkl,38,NL.1RA,12
1,112,iLc N1,Z,rlO,NL,lrIl,121,142,It+
2
03 05
Nu02,Nk,NL.lnnI.ILi.,Il.s,12+ NOZ,Ni,NL,141.ILA,Irr(,IL4 N82,I1o,r,L,1%1,ILi,1nC2.12c
3
03
~
os
NRZ.N4,NL,I4l,1ui,1R4,IZ~
63
63
M82,HnI,.,t,icu.I2n.,IA2,l22 082,38,52,Iii.iZl.I*2,IiA NR2,9t,41,18l,11L,102,AZl NRZ,Iir(,n42,l1,L,ILL,Ii14,1u2 NI1l,N8,NL,111,LLI,Iil.112 NRZ,N*n,NZ.I,nI.IL1,1’2,1Zc
05
+3
o7 53 SI
4*3 uS
ru8Z,kA,N2,inrl,I2i,I+2.IZ.~
7’,
?u
N02.44,r12.14i.12i,154,121 N8L,MR,1,44.I1L,12n.,ilo2,122
7+ /3 i~
P,82,Nr,nMZ,Ioi,I21.1,n2,iZs NR2,N
is
5,NL,1k1,ILi,162,Ii.. NRZ,N*,91,1,*1,12i,I*4,122 482,N6,N2,1*l,I22,IIrd,I22 NRZ,N8,niL,III,IZn.,152,I14
NkL.r48,IIL,lli.1EA,102.112 N82,Nc,51,I,IA,I24,I’r2,12+ NR2,98,52,101.ILi,i42.1L4 nlR2,N0,n0L,irl,1L~.,t42,iL1
—.ou9111—01 —.173b61.uU —.144791,1,1 —.184481*10 —.171201+40 —.154941*00 —.640241—01 4
0:,
30
—.o73Lu+—~i .jsls.+—U+ 7+000 —.L3nn3,1*sn —.+no4n.41*su —.i/32 —.+449)t*u3 —.5+USoi,3i mu
47
0:, 0+
s:,
65 5, 53
33
+4
i
33 33 33 uJ
64 6’, 04 i4
2
33 33 33 ,3 33 ii 33 :,3 33 33 33
ii uI 30 2’I 4* 27
31
+0 +4
15 43 42
55 SS 05 sin
3 4 3
3c 32 32 ii 3i. 34 3~ 32 32
4 3
45 14 43 24
04
33 33
2
*1
1
01
4
7
III
03
33
0
i 1
i23 124 is,
+
12/
03
33 33 33 33 33
3
048
4) o5 03 55
2 2 2 4 4
73 —.354,11—04 0i.1u—ui —. +7os*.+,.uO —.l3:,721’sO —.102035*14 —. L4445c*30 —. 303932—.
3
33 33 33 33 33 .13 ii .13 uJ 33 33
1+u +41
—.
2
.12 ol 35 49 2+
1
no 7 6 5
4
£
4
3
4
3
i
4
+
2
1
1
1434+1—+i —.31353+—Ui .lT543lOUnn —.+40U5u*lu —.13125*054 —.+3*97+*niu —.235314—ui —.
369
i
03
17 25
stability code
—.16495+—Ui —.37903+—li —.140s01000 —.14)no2+*oO —.i79s9+oUs —.127o4+*s~ .O13?ns11c
3ZbilInnI —.254241—41 —.944604—UI —.162111*50 —.I5o:,3+*00 —.176421000 —.110441.00 .41U151—02
—. —.104330+01 —.18355t+uU —.lSSnnUEOuU —.A73431*Uu —.iUOi5E*Ul..
—.391731—01
—.468811—Ui
—.147051.v0 —.184541*00
—.
—.159221.00 —.169541*00 —.990651—01
Ilji3louu
—.185051*00 —.l63i6+*uU —.165221.0+ —.68481+—ui.
—.53930101.
—.ii9USE.00
—.185188*04 —.150831*011 —.160371*04 —.77u191—Ul
1+
FIT u.
3.54o+—u3
,.i57,.—UI
u.
ni. 4.1331—is
4.
0.
4.
0.
5.7101—4+ 6.450—si 7.9431—ui 9.1.151—Ui 1.o251*uu
5.83I~—S1 5.98.11—ui 5.04s1—ui l.2~nis—4L 1.0301*14 1.14*1+03 1.26,+.1.u
in.io4+—uL 7.077+—SI 0.lu+S1 9.34 c15 1 1.J47n.’ss 1.16150+0 +.274n.~1u
~.A47t—nii 7. Alit—ui no.32u1,u1 1. 4551—uI ~~1,591*o5 i.1721.so i.2uSt,00
0.Li++—Sn. 7.3431—UI 3.4,1,2—ui 3.10 ++~0n. i.ulU+*sJ 1.1331*00 1.4.7+0.5
s.204+—ul. 7.4061—Ui no.0.91—ui 9.5911—01 A.10i1*uO 1.li34.*sU
1.
4.333+si
4.Lb/+—ui
u.
0. 4,1331—lu
4.101+—niL u.
.n. 4.1572—si
3.1*44*—ni u.
.,.
4.,o71—0i
5.4311—si o.9831—ui. ,.u’1~i—s1 1.4235—ui i.53a1*nnl 1.149*055 1.2531*40
3.9441—uI
o. US7+~niI 7.1954—si 1.32~t—u1 9.4531—ink I.+:,9100s L.u?4s*Un.n 1.2451*oU
—1.8521—01
uc7L—uI
7u01—11
4.1841—01
2.4881—lU
0.
6.397+—Ui. 7. S3Uc—Uj 6.602+0I 9.795+—UI A.09310UU 1.2Uoi*UU
0.5101—UI 7.6431—Ui no.7754—01 9.9081—UI i.i041~0U n..2171*UU
6.6241—61 7.7561—UI 8.8891—01 1.0U21+uO i.i151*0J 1.2291.00
0.7371—UI 7.8691—01 9.0021—Ui I. 013E.UU 1.127E*uO i.24u100u
0. U.
—2.72,1—11
3.9811—11*
2.4891—l0
In.
u. 0.
—1. 7bU+—11
4.184*—UI
2.5361—10
u.
6.3971—LI 7.5341—UI *.6o21—01 9.7951.—li 1.0931*00 1.4151*1n4
0.5151—Ui ?.5431—uI ,.77~1—ul 9.9041—41 i.1041*uU 1.2171.00
6.6241—Ui 7.7561—~u1 6.8841—UI 1.0021*00 1.iI5l,UU I.2291*UU
6.7378—01 /.Ut,9101 9.0021—UI 1.0131*1*0 A.1271+uU 1.24U..Un3
3.
—2.
108
1..134tnosnn l.cSit’Us
AX1S~
39
FIT 4.157+—un. 0.
FIT —1.0541—Li 0.
3. 98
lUsU
0.
104
XIS~
5. 7lnot—Ui 6.OSuE—us l.9031—ui 9.1151—01 1.u231*uu 1.I38i*IU i. 2,11 *40 39
7.s71+—si. 1.4091—ui 1.344+—nil l.u47+*uu i.I5+1*niU 1. 274+ *sn3
silt—ui 7.3s3+—3i .3.43b1.n3n. 3.500+—nil i.07n12*i. L.1,~~0nU i • 2* l+*us +.
5.4*141—Ui 7.4101Ui no.349.—Ui ‘3.63+L—Ai 1.0011*00 1.1*51000
370
R. Gruber etaL
/ ERA TO stability code
JP~I
01*
OsIIT*
.252970~2 .s19i1943 .43235833 .323b356. .44859254 .4,224729 .36775902
.s777uo 7u
.1.1933887 .11*770044 .4110+534 .1)423914 .u159s7,.3
.01619807
• A~, A3inon.,5 .14400414 .+1924353 ..1397u97, .0010+445 • 1~0S35 ,o7 .~,U4ol+4 .0360145.. .s:,5js3:n+
.34244139
.32434 312
.258358~7 .07695517 .13904592 .01,741350 .34227044 .u77064ld .n.ni89.,177
0,14*5,0 .30448952 .1i5,314~ 1.43404410 +.*n7+Ijodn, i.,Ssl+.+, 3.nn3u3493s 2. ni9i43’,3:, i.9U3o4o’~5 1.24*/483+ .0in)~,,+ii .015344+0 • 4bL+531u .3o7,345no
.31235042 .49779144 .1 1743340
*70 .23403433 .2469*1431 .4,045223 ..1391*Dull 44nL97b00 .i49.41,n, • 1337*347 .1l1i5323 .2450)7+4 .44473+13
.s1734353 • oui’i SFn.o .01)09310 .5,3357104 .00290)40 .
04243171
3~S1 011*10 .35497573
.i30,..3u,. • 1447iu 15 • 440695 As
1.55341343 1.s313 *144 4.13734273
.94985319
.3117/3*2+
2.5)5L,,207
.036*444n.3
.27*14,40)
.05494+31 .351734951
.2743s..4, .14’llt,nns .5.17,/,43
3.3194*337 4.3, liiuic 3.337.7151 +...JuUt.11i
.1,54)8,1
.212o os41
..5 4nn~3u,~o
.n.15,000’0
...1o7’,191
1.13122.7+
.6)911744 .03099310 .u14u4s9S .uni4550o4 .u3i47 192
.i,79+SSo .+,+53/37 .i42ci,5, .43773+14 • 35+4443
.17j4+oz5
.0)104101
.5,4343.7
•351212~5
+0070333
.14234117
.133+473’, .++i70o/+ .44+79.04 .40171531
431110
.311351)7 .3334325’i .46 333i05.
.b018477i .35703045 .75903984
:,)94o
72
.35614/94 .21612533
.1243.301 .ni9I.0usI3 .v229YU*4
22 53io37 .09000+10 . :,r,64,’, *1 .3132,3/7 .llSuO4snl •
.4,nilnoin’1’, • 3177/37+ .3,11+0 + 3 • 33323014 .3)277018
.i)639*
:,u
*00
—.+n77,8751, 1303471
—.53393263
.s3403+uO 4j ,344 .s~i75)+4 .4110794’, .o..’n,n*.45i
—.75s4641*l —i.-1
i.24354300
—...
—ILn.Ln.n33334
—.33u,2:,2~ —.A5i*ils
.51)54070
—..,in.,4711n4 —.4iu74)/5
.un..?373o . ,~I 3,375
—.5471?i41 —.01.11,511
.s24n,73n—.04141240—. =
—.039307/3 —1.12795931 — 1.7,409908
—3.35069932 —5.307541*4 —4.o.b77~+I .03430347 1.11777126 .9432459k .05)24958 .8937s99~ .si9iu:,Sl .oI /o36o1 10307)70
.2
47+
4.01,52511’
.44+365*43 .3,533001
~i.lcs~4S4*.
.7421n339)
on’s
•
1,3+795
—.4344915, —.llIs+550 —.4)7,44c0
+1
s.+4o3741,5
4.
—.1,5333,57
—4.419753,’, —4.44271413 — 1.4*l+n+4j.~ .301*04277 .02533033 .7078)53+ .o2oio~,O • 4952093*1 .43036u21 .4)040794 .33+936u3
—2.01104157 3,147124 —4. 5007 i27 /..ni~711o9 —sO.31349308
4;’.j/.~~4
1. .oo4i/5/ 3.+3325i~4i 1.73301351 4.15,34,05 3.74 170590 4.42130 30,
8.34+.)no56 —5.iSsH13 —d.O4ouul3ln —3. .,so,94352 —2.us2o?1*53 —.77niuo93:, —.34.2,33*4 —. 17413)13
—.39886044 —1.40227441 —Z.j39159no6
.2305491, • 2s4.4737 • ci’,9143s .47043331 .,.20n.n.4u9 .u,720o4o .iuliUuq+ .,0
3/Si 019
*0*
—l.*1lo4n,HiO —4.24,37233 —3.1+7,3750 77
07*
~01540,
• 3205514
*00
l,uss’*/
.o
41314+0 .42450334 • 11444905 .s127 4S9i .u*423i .,9* • ~,21,418+ .ui2:,*7n,3 .45319134 .50494023
.01213205
.43541050 .2*3,suo7 .1372,0,9
.593,3174 • 74n,59+71 ./411+4,2
.004445+2
.335+0543
.70*3933+
.,u052343 .u405s4+u
.0o715211
.u0319775
.44+99+35
.s2i29602
.uu47nIsio—.
444
3.L
00*
ooun7inl
—.,.l,Su,h, ).977
191*31
—9.35045111
—43.bilOuiSO —io./6~7u764 —4.3+11*4,51 —1. 7004.23’u .00/0)+43 —.2354,’,91
—I. uo4ISo*49 —o.72030i11, —2.43*87055
—4.51174110 —7.j3so2lAA —4.49i33079 —. ~5
949434
2.44.35037 £ .o2204u 70
1.40344956
—.14303933 —.1,043)7)7
U5+3~36s —.0394114) —.03750+37
—.73131183
1.11438241 i.uU17ll144 .3,2831)5
.31411291 —.i737454u
R. Gruber et aL / ERA TO .3/SI *01070 .78370.43 .50340419 .64240694 .93524886 1.1955059b 1.01014193 .774u8i6o .498595,3
.3371.4+4 .42714344 .49054400 .,121.8943 .5I52’18u2 .5425434*3 .47490925 .43704463
1.19322413 l.u7inA*349*
.44481407
—14.4333957*)
3.443,?slb
4.l4,~oO3u4 3.23+24574 2.133,0947
.,.I9+077’, .1444443/ 759/ .nn33. .sno),l5l
—10.09991129 ~0./54)7832 —4.33881363
in.S)535044 i.lbioSmju .u5s3533, .04231542 .7o,~7o5i
.u24213o3 .0117501,4 .uu354:,,o .1u59?u?3 .ouu14~s9
.05071,11
.09453965
.017547u9
.40,o1733
.n,3.44355
.3321.377 • 47944*3, .5255025, .34533004 .43997745 .4,12.31,9 .40572,35
.141557*30
.40143275
.3370,063
• iS7sosY3 .2U5+,~~5 .15759,11, .12394313 .14253424 .u300304*, .o37niuS64
4IIX+1*4
—1.30129356 —.))5u5335 —.22471347 —.u9d15sU~ —.44072525 —.424u.408
—2.55153044 —3.49055+63
1.6560117? 1.74385126 ..32177047 1.31189787 1.1320+786 1.03615253
•03I~0jo/ .514o21.,ni—.
—.30450143 41103347
.9703911+ .94697030
.s1327+55
—.04979707
—.23622044
3751
1.0*1110913 .71946347 .7947444o 1.236nn1543 1.659059*3, i.529145uu 1.12502808 .73740044 .49059839 .3542Uc78 .28231034 .4394459? .41294985 .19819342 .19122520
—5.9755’8431 —0.91630780 —4.39U2999, .27487002
3.04979543
.23441101 .43919451
*00 —1.551.31163
2.44970044
.u72159o7 .06816438
*~)17T*
04*
3.973,7339 —4.73391330 —7.04934940
.1135+474 .175/4*4+ .47/3*7,3
.33357544 •3sn.,’9129+ .3149004, .4312908s .27i+7~o4
371
1.i
*70
.30871534 .20211,101 .141900)1 .10747535 .04534318
*1*
=
stability code
•
9*’4
1.47074414 .10354301 ~s4212z34 3.3383455’, 4.4,’,2544c 3.54317679 2.o0~1541 0.3742171) • ,71s 718+ .3/74)85+ .397.379*1’, .2+447923 .21I:,2nio0 • 173 + 7+nii .11227539
.u13741sn. • o45*3~L4 .160)4,73 .i9177023 .I343s993 •i35o2705 .o~no7005., .40204045 •uc4bOUis .014,4004 .u07U43oA .0543+724 .112,3430 • uu1b3399 .0ui365,4
.2=74?s04
.uL~55o54
*40 —4.07+87444 —5.8~7~2420 —9.~7/+o8d1 —15.89533192 —n.no.7715,135 —.u.1nnU39523 —3.93857077
—1.35/5in397 —.51340314 —.40.73023 —.0,027047
—.u3?o9913 —.51504552 —.337030.9 —.uo4sSSOo —.o)541595
*48 —2.31969036 —3.66089913
—~.5)420n373 —20.40016102 —11.97)12,73 —4.9550310o —.48092042 A.134u0229 1.15106791 .94050834 .~79lbI52 .011431,45 .51349331 .4,679041. .93139576 —.44,o*31i2