ORCO: A numerical tool to study the radial diffusion of runaway electrons in tokamaks

Computer Physics Communications 156 (2003) 95–107 www.elsevier.com/locate/cpc

ORCO: A numerical tool to study the radial diffusion of runaway electrons in tokamaks R. Sánchez a,∗ , J.R. Martín-Solís a , B. Esposito b a Departamento de Física, Universidad Carlos III, Avda. de la Universidad 30, Leganés, 28911 Madrid, Spain b Associazione Euratom-ENEA sulla Fusione, C.R. Frascati, C.P. 65, I-00044 Frascati (Rome), Italy

Received 31 January 2003; accepted 30 July 2003

Abstract ORCO is a new code that estimates the spatial structure of the radial diffusion coefficient for runaway electrons in tokamak geometry. In real experiments, the location of these electrons can be detected by measuring the time evolution of their fast electron bremsstrahlung (FEB) emissivities, usually integrated along several lines of view. ORCO uses a Levenberg–Marquardt algorithm to adjust the free parameters of a generalized transport model to best reproduce the time evolution of these temporal traces. A possible future application for this type of calculations is to use them as indirect probes to test theoretical models of turbulent transport driven by stochastic magnetic fields in tokamaks.  2003 Elsevier B.V. All rights reserved. PACS: 02.60.Lj; 52.35.Py; 52.55.Mc; 52.55.Kj Keywords: Runaway electrons; Tokamaks; Magnetic turbulence; Diffusive transport; Non-orthogonal coordinates; Levenberg–Marquart optimization

1. Introduction In tokamak discharges, runaway electrons can be generated whenever a large toroidal electric field is induced. This is the case, for instance, during the ramp-up or the quench of the toroidal current. This electric field continuously accelerates electrons with initial energies exceeding some critical threshold, that is set by the value beyond which the collisional drag can no longer dissipate all the energy that electrons gain in the electric field. These electrons are known as * Corresponding author.

E-mail address: [email protected] (R. Sánchez). 0010-4655/$ – see front matter  2003 Elsevier B.V. All rights reserved. doi:10.1016/S0010-4655(03)00439-9

first generation or Dreicer runaways [1] and would in principle accelerate with no bound. However, a second energy sink associated to the increasing radiation losses, that eventually dominates at higher velocities, sets the final runaway energy [2]. Once a sufficiently large population of these Dreicer runaways builds up, secondary generation or avalanching runaways can be produced as well as a result of the collisions between Dreicer runaways and thermal electrons [3]. A common feature to both types of runaways is that, since they scarcely collide with thermal particles, they diffuse out of the tokamak by following the magnetic field lines. (This is in contrast to what happens for lower-energy electrons [4].) Thus, the quality of their confinement could be interpreted as an indi-

96

R. Sánchez et al. / Computer Physics Communications 156 (2003) 95–107

rect estimation of the levels of magnetic turbulence and stochasticity present inside the tokamak. Indeed, it has been shown that, in the extreme limit of a fully stochastic field in toroidal geometry, the runaway diffusion coefficient can be estimated as [5,6]: D(r, v ) = Dst v ,


|bmn (r)|2  m −n , δ Dst (r) = πR q(r) B02 m,n

(1)

where v is the runaway velocity, bmn is the amplitude of the resonant field mode located at the minor radius r∗ that verifies ι(r∗ ) = m/n (ι = 1/q(r) is the rotational transform of the tokamak-like configuration) and B0 is the non-perturbed toroidal magnetic field. Thus, some knowledge of the v -averaged value of D could be sufficient to estimate the levels of magnetic turbulence present of the device. In the opposite limit, for a magnetic toroidal configuration with perfect closed magnetic surfaces, the diffusive coefficient would be given by the neoclassical one [7] (evaluated in the collisionality and velocity conditions relevant for runaways). In a real tokamak, regions of stochasticity would probably alternate with good-quality closed surfaces, and the real diffusive coefficient should lie somewhere in between these limits [8]. (More recent discussions on these important issues can also be found in Refs. [9,10].) In the last fifteen years, it has been possible to obtain estimations of an spatially-averaged, velocityaveraged radial diffusion coefficient for runaways, from a variety of diagnostic systems, in several tokamaks [11,12]. These results have been then used to estimate the levels of magnetic turbulence existing in these devices, as well as to test to some extent various theoretical models of stochastic transport [13]. However, no information could be gained about the spatial distribution of magnetic turbulence due to the radial averaging. In this paper, a code is presented that attempts to attack the spatial deficiencies of previous approaches in a numerically efficient way (some comments will be made later on the remaining issue of the velocity dependence of the runaway diffusivities). For this purpose, the ORCO code (an acronym for nonOrthogonal Runaway Confinement Optimizer) implements a generalized transport model that tries to capture the basic physics of runaway production and

transport in tokamaks while keeping information of their spatial structure. The model considers both diffusive and convective transport as possible processes available to drive runaways out of the tokamak. Information about the spatial dependence and relative importance of these transport channels is stored in a set of free parameters. The values of these parameters that best match the experimental runaway observations are then searched for numerically by means of a Levenberg–Marquardt optimizing loop. In the version of ORCO presented in this paper, the experimental data to be matched consist of line-integrated measurements of the fast electron bremsstrahlung (FEB) radiation emitted by runaways during their scarce collisions with other particles. This radiation can be routinely detected in several tokamaks such as the Joint European Torus (JET) [14] or the Frascati Tokamak Upgrade (FTU) [15], and related to the runaway electron local density by [11]: IFEB (r , t)  Cr nr (r , t)ne (r , t)Zeff (r , t),

(2)

where nr denotes the runaway local density, ne is the local thermal electron density, Zeff is the effective charge number and Cr is an unknown constant. However, it is important to keep in mind that usage of ORCO need not be restricted to the geometry of this particular type of diagnostics, thanks to its highly modular structure. In order to minimize the errors that are unavoidably added by the numerical optimization to the already large experimental uncertainties, several issues have been taken into account in the transport model within ORCO. A first one has to do with the motion of the plasma column relative to the fixed lab frame in which the FEB data are measured. This motion may cause that distinct spatial locations contribute to a single temporal trace of FEB data. A second issue has to do with the deformation of the plasma equilibrium geometry, that can change as a result of the variation of the plasma currents while data are being recorded. ORCO tries to minimize these problems by using to its advantage other available experimental information regarding the position and shape of the plasma. This information is accommodated within the transport model by means of a time-dependent coordinate system, that constitutes one of ORCO’s most characteristic features.

R. Sánchez et al. / Computer Physics Communications 156 (2003) 95–107

The paper is thus organized as follows. First, the general structure of ORCO is described in detail across several sections. In Section 2, ORCO’s internal timedependent coordinate system is introduced. Then, in Section 3, the FEB diagnostic system existing in both JET and FTU, to which the transport model is currently coupled, is described. In the same section, the intersection of the lines of view of this system with the plasma column are analytically calculated. Next, ORCO’s transport model equation is derived in Section 4. In Section 5, the numerical scheme to integrate this model is presented. Finally, in Section 6, some details on the Levenberg–Marquardt optimizing loop are given. After the presentation of ORCO, its application to a test case (a limiter discharge from the JET tokamak) is described as an example in Section 7. Finally, some conclusions will be drawn in Section 8.

97

2. ORCO’s geometry and coordinates As it was mentioned in the introduction, it is important to minimize any error that the numerical estimation of the runaway transport model might add to the uncertainties already present in the experimental data. A first problem that needs to be dealt with is the motion and the changes in shape of the confining magnetic structure during the discharge. As a result of this motion and changes, the regions of the plasma intersected by the fixed (in the lab frame) FEB diagnostic lines of view (see Fig. 1, that shows a sketch of the FEB diagnostic at JET), along which the FEB emissivity is integrated, do change with time. Therefore, these shifts can easily interfere with the correct interpretation of the experimental FEB data if not properly taken into account. As a result, the reliability of the runaway diffusion coefficient

Fig. 1. Fast Electron Bremsstrahlung (FEB) radiation diagnostic system at JET. Horizontal channels are numbered #1–10, and vertical ones #11–19.

98

R. Sánchez et al. / Computer Physics Communications 156 (2003) 95–107

obtained by ORCO and the conclusions derived from it might be seriously compromised. ORCO faces this problem by means of a nonorthogonal time-dependent coordinate system that allows to use to its advantage all the available experimental knowledge of these changes during the discharge (see Ref. [16] for a thorough introduction to generalized coordinates). In it, each point in the magnetic configuration is described by three coordinates, (F, φ, θ ), where F labels a magnetic surface (going from 0, at the magnetic axis, to 1, at the last closed surface (from now on, referred to as LCS)) and θ and φ are two poloidal and toroidal geometrical angles (varying between 0 and 2π ) that locate the point on each magnetic surface. A mapping to these coordinates can be easily constructed from the usual cylindrical coordinates (R, φ, Z) if it is assumed that the LCS can be well described by an ellipse of minor and major axis respectively given by a and b (see Fig. 2). To do it, a first functional F ∗ is defined as follows (this function is a slightly modified version of a parametric solution of the Grad–Shafranov equations for a high-β tokamak geometry [17]): 

(R − R0 )2 (Z − Z0 )2 F (R, Z) = + −1 a2 b2   (R − R0 ) + 1. × 1 + k2 a ∗



This functional verifies that F ∗ (R, Z) = 1 corresponds to the elliptical shape assumed for the LCS, with its center located at (R0 , Z0 ), as long as |k2 | < 1. The remaining undefined parameter, k2 , is then prescribed to match the location of the magnetic axis. Indeed, it is straightforward to locate the minimum of F ∗ (R, Z) at (Raxis , Zaxis ) given by:  a( 1 + 3k22 − 1) Raxis = R0 + , (4) 3k2 Zaxis = Z0 . Defining now a normalized axis horizontal shift respect to the center of the LCS,  ≡ (Raxis − R0 )/a, and inverting Eq. (4), it is easily obtained that: 2 . (5) 1 − 3 2 The constraint |k2 | < 1 thus limits its applicability to discharges with a maximum Shafranov shift ||
1/3. This condition is however fulfilled by most [moderate-β] limiter discharges. Finally, after normalizing F ∗ appropriately, a suitable expression for F (R, Z) reaching its minimum and vanishing at the magnetic axis is obtained: k2 =

F (R, Z) ≡ (3)

∗ F ∗ (R, Z) − Faxis , ∗ 1 − Faxis

where the value of F ∗ at the axis is given by: ∗ Faxis ≡ F ∗ (Raxis , Zaxis)    2 + 1 = − 2 . 1 − 3 2

Regarding the θ coordinate, it is defined by:   Z − Zaxis θ = tan−1 , R − Raxis

Fig. 2. Sketch explaining the relationship between cylindrical (R, φ, Z) coordinates and ORCO’s (F, φ, θ ) coordinates.  is the magnetic axis (MA) shift with respect to the last closed magnetic surface (LCS) center.

(6)

(7)

(8)

that corresponds to the geometrical poloidal angle centered at the magnetic axis. Finally, the toroidal angle φ is trivially defined to be the same that the azimuthal angle of cylindrical coordinates (see Fig. 2). Since they will be needed in later sections, we will also give here the prescription to obtain several quantities related to the metric of the new coordinate system. √ In particular, the Jacobian of the transformation, g:   ∂F ∂θ −1 ∂F ∂θ √ − (9) g=R ∂R ∂Z ∂Z ∂R

R. Sánchez et al. / Computer Physics Communications 156 (2003) 95–107

and the contravariant metric element g F F :     ∂F 2 ∂F 2 FF + . g = ∇F · ∇F = ∂R ∂Z

(10)

In summary, the shape and evolution of the equilibrium geometry is totally determined in ORCO’s (F, φ, θ ) coordinates by five functions of time: a(t) and b(t) that give the change of shape of the LCS, R0 (t) and Z0 (t) that give the motion of the center of the LCS, and (t) that contains any possible change of the Shafranov shift. All of them can be routinely obtained from the reconstruction of the discharge under analysis with the EFIT code [18]. [EFIT is a MHD equilibrium solver that, thanks to the addition of appropriate response functions for existent diagnostics (such as magnetic probes, poloidal flux loops or Motional Stark Effect measurements), is able of reconstructing the magnetic geometry of a given discharge fast and efficiently.] It is also relevant to say that the representation presently chosen for F is only valid for ellipsoidal plasmas, as those encountered in limiter tokamak discharges. The reason for this choice is that, up to now, the code has been used to analyze limiter discharges from the JET tokamak and, in the near future, it is also to be applied to FEB data from the FTU tokamak. However, other representations could be easily implemented to parametrize different geometries (for instance, D-shaped divertor plasmas). As a matter of fact, envisioning the future use of the EFIT equilibrium solution (or that of any other MHD equilibrium solver) within ORCO would not pose a particularly difficult problem.

3. Description of the FEB diagnostic system The FEB experimental diagnostic available at JET [19] consists of two cameras that provide line-integrated temporal signals of the X-ray FEB emissivities along ten vertical chords and nine horizontal chords that intersect with the plasma at constant toroidal angle (see Fig. 1). It collects data in four different energy windows (in KeV): [133, 200], [200, 266], [266, 333] and [333, 400] (see Ref. [11] for details). A similar FEB system (on loan from CEA-Cadarache) consisting of one camera with 17 chords was also used in FTU during the period 1998–1999 (see Fig. 2 in Ref. [15]). Currently, a new FEB system consisting of

99

two cameras (vertical and horizontal) is being installed in FTU and will be operative in the near future. In ORCO, all the FEB chords are parametrized by straight lines using the usual cylindrical coordinates: Z = mk R + ck ,

(11)

so that a pair of values (mk , ck ) represents each of them. Then, the intersection of the kth chord with the changing plasma is calculated analytically by looking for solutions of the following system of equations:      Z − Z0 (t) − mk R − R0 (t) = Ck (t)  2 2 (12) [R − R0 (t)] [Z − Z0 (t)]   + =1 2 2 a(t) b(t) where Ck (t) ≡ ck − Z0 (t) + mk R0 (t). A necessary and sufficient condition for the actual intersection of the kth chord with the plasma at time t can then be stated as: k (t)2 ≡ b(t)2 + m2k a(t)2 − Ck (t)2  0,

(13)

and the coordinates of the two points of intersection of the chord and the LCS (when Eq. (13) is satisfied), k k k k (R− , Z− ) and (R+ , Z+ ), by: k (t) = R0 (t) − R±

and k Z± (t) = Z0 (t) +

mk a(t)Ck (t) ∓ b(t)k (t) b(t)2 + m2k a(t)2

(14)

  k 2 b(t) a(t)2 − R± − R0 (t) . (15) a(t)

k k , Z± ) will prove very The analytic knowledge of (R± useful later, since ORCO has to compute at each time the line-integrated FEB intensities from the solution of the transport model to compare them with the experimental FEB data.

4. Runaway transport model In ORCO, a single partial differential equation is used to describe the radial transport of runaways out of the tokamak. It includes a runaway source term as well as possible diffusive and convective loss channels. However, an extra term must also be included due to the built-in time dependence of ORCO’s coordinates through a(t), b(t), (t), R0 (t) and Z0 (t). This dependence makes the local volume differential, dV =

100

R. Sánchez et al. / Computer Physics Communications 156 (2003) 95–107

√

g dF dφ dθ , to depend also on time. Therefore, the local runaway density may change with time in these coordinates, even in the absence of local sources and/or sinks! However, the fact that the local number of runaways must be conserved implies that √ nr (F, θ, φ, t) g(F, θ, φ, t) √ = nr (F, θ, φ, t + t) g(F, θ, φ, t + t). (16) Using this equation, the extra term needed can be easily estimated by calculating the temporal change of runaway density in the absence of any source or net flux of runaways: nr (F, θ, φ, t + t) − nr (F, θ, φ, t) ∂nr = lim t →0 ∂t t √ nr ∂ g = −√ (17) . g ∂t Therefore, ORCO’s runaway transport model will be derived from the following three-dimensional (3D) partial differential equation: √ nr ∂ g ∂nr +√ ∂t g ∂t  = Vc · ∇nr + ∇ · Γrd + Sr (r , nr , t), (18) where Γrd is a diffusive runaway electron flux and Vc is a runaway convective or drift velocity. Regarding the runaway source, Sr , it includes the two generation mechanisms mentioned in the introduction: nr p . Sr (r , nr , t) = Sr (r , t) + (19) τs (r , t)

code widely used within the tokamak community that works on a two-step basis: first, profiles are advanced according to a prescribed transport model; then, for fixed profiles, a new equilibrium solution is computed. By using appropriate initial conditions provided by experimental data, many quantities of interest can be easily computed. Information on the most recent version can be found in Ref. [25].] Regarding the diffusive runaway flux, it is given by Fick’s law:   ∂nr ∂nr ∂nr Γrd = D(F, φ, θ ) ∇F + ∇φ + ∇θ , ∂F ∂φ ∂θ (20) and its divergence is equal to:     1 ∂nr ∂nr ∂ √ ∇ · Γrd = √ + gF θ gD g F F g ∂F ∂F ∂θ    ∂nr ∂nr ∂ √ + g θθ gD g F θ + ∂θ ∂F ∂θ   ∂ ∂nr √ + (21) g φφ gD , ∂φ ∂φ where g ij ≡ ∇i · ∇j (i, j = F , θ or φ) stands for the different upper metric elements of the coordinate transformation. Since we will restrict ourselves to studying the radial transport of runaways, Eq. (18) is reduced to a one-dimensional (1D) equation by surface-averaging. The surface-average of any arbitrary function A(F, φ, θ ) is defined to be:

p

Sr gives the Dreicer generation rate. Its original definition can be found in Ref. [1], and a better analytical estimate in Ref. [20]. However, the one used in the code is based on the calculations reported in Ref. [21], that includes relativistic corrections. On the other hand, 1/τs represents the characteristic production rate via avalanching, first estimated in Ref. [3]. The expression actually used in ORCO is that of Ref. [22], which is based on the calculations first carried out in Ref. [23]. Both of them are functions of available experimental data such as the electron density, ne , electron temperature, Te , effective charge number Zeff and electric field, E . [The electric field is obtained either from the experimental loop voltage or by reconstruction with the JETTO code [24]. JETTO is a transport

1 A ≡ 4π 2 G

2π 2π dφ dθ

√

gA(F, φ, θ ),

(22)

0 0

where the function G is defined by:

G(F, t) ≡

1 4π 2

2π 2π dφ dθ

√ g,

(23)

0 0

that can be interpreted by realizing that 4π 2 G(F, t)F is the volume enclosed between the magnetic surfaces labeled by F and F + F . The surface-average of Eq. (18) then yields:

R. Sánchez et al. / Computer Physics Communications 156 (2003) 95–107

 √  nr ∂ g ∂nr  + √ ∂t g ∂t     1 ∂ F F ∂nr F θ ∂nr +g = G D g G ∂F ∂F ∂θ  + Vc · ∇nr  + Sr .

101

meshes are built using temporal and spatial step sizes respectively given by: t ≡ (24)

Then, the runaway density is rewritten as nr (F, φ, θ ) = tnr  + nr,1 (F, φ, θ ) and it is assumed that |nr,1 |  nr . An effective surface-averaged radial diffusion coefficient is defined as well:

tf − t0 Nt − 1

and F ≡

1 , N −1

(29)

where tf and t0 are the final and initial times; Nt and N are the total number of points respectively in the temporal and (full) spatial meshes:

(25)

tj = t0 + j t, j = 0, 1, . . . , Nt − 1, k = 0, 1, . . . , N − 1, (30) Fk = kF, Fl+1/2 = Fl + (F /2), l = 0, 1, . . . , N − 2.

This allows to reduce the transport equation, after several straightforward manipulations, to its final form:

The solution at time tj is represented by the vector:   j j j j Nr ≡ (nr )1 , (nr )2 , . . . , (nr )N−1 . (31)

D(F ) ≡

g F F D . g F F 

∂nr  nr  ∂G + ∂t G ∂t  ∂nr  1 ∂ Gg F F D = G ∂F ∂F  p n  ∂n r r + VcF  + + Sr , ∂F τs 

(26)

where the surface-average of the F -contravariant component of the drift velocity is defined as:  F Vc ≡ Vr · ∇F , (27) and that must satisfy that:   lim VcF = 0,

F →0

(28)

since ∇F vanishes at the magnetic axis. When examining Eq. (26), it is clear that the complete determination of the transport model only requires the finding of two free functions of the surface label F : a surface-averaged convective velocity perpendicular to the magnetic surfaces and an effective surface-averaged radial diffusion coefficient. The optimizing loop inside ORCO will search for them by requesting Eq. (26) to be able to reproduce the experimental FEB measurements (see Section 6).

5. Numerical scheme To integrate Eq. (26) numerically, a stable and accurate iterative algorithm has been derived by using a second order Crank–Nicholson differencing scheme [26]. The temporal and spatial (both full and half)

Note that the values of nr at F0 and FN are fixed for j all time by the boundary conditions chosen: (nr )0 = j j (nr )1 , ∀j and (nr )N = 0, ∀j , so they do not need be included in the iterating scheme. The resulting iterating scheme can be set in matricial form: j +1

Aj +1 Nr

= B j Nr + C j,j +1 , j

(32)

where A is a (N − 2) × (N − 2) tridiagonal nonsymmetric matrix given by:  j −1 3 Gk t j Akl = δkl − − j j 2 2Gk 2τs k  t j + g F F k−1/2 Dk−1/2 j 2(x)2Gk  FF j + g k+1/2 Dk+1/2  FF j g k+1/2 Dk+1/2 t −δk,l+1 2 j 2(x) Gk  F − xVc k  FF j g k−1/2 Dk−1/2 t −δk,l−1 j 2 2(x) Gk  F + xVc k ,

(33)

being δj k the usual Kronecker delta. B is a second (N − 2) × (N − 2) tridiagonal non-symmetric matrix related to A through (I is the identity matrix):

102

j Bkl

R. Sánchez et al. / Computer Physics Communications 156 (2003) 95–107

  j +1  j −1  1 Gk + Gk j = 3− δkl − Akl . j 2 G

(34)

k

Finally, the independent (N − 2)-vector C j,j +1 is given by: j,j +1

Ck

=

t  j 2Gk

j +1

Gk

j +1

S p k

j j + Gk S p k .

(35)

Notice that in all the previous definitions, integer [semi-integer] subscripts denote evaluation on the full [half] spatial grid, while superscripts refer to the temporal mesh. The Ax = B linear system of Eq. (32) can now be easily solved by using a standard LU decomposition routine for tridiagonal matrices [26], since j +1 C j,j +1 is independent of Nr , being completely determined at all times from the experimental data. Note j as well that all Gk are known analytically at all times. 6. Levenberg–Maquardt FEB scheme The optimization module is the part of ORCO that explores the free parameter space in search for the D and VcF  functions that best reconstruct the experimental FEB measurements. This search is carried out using an optimizing algorithm based on a Levenberg–Marquardt loop [27]. The way it is done follows a scheme that has already been proved very successful in the design and optimization of compact stellarators [28–30]. It goes as follows: first, a target function is built that the optimizer will try to minimize with respect to a vector of free parameters, p = (p1 , p2 , . . . , pMp ). This vector includes both D and VcF  (using a representation that will be discussed below) as well as the unknown constant Cr appearing in Eq. (2). The target function mainly carries information about how well the experimental data can be reconstructed from the model. But some other constraints that must be satisfied during the optimization are also included (for instance, the positiveness of the diffusion coefficient on all surfaces). Thus, it is built as a sum of Mt partial target functions χk : χ 2 (p)  =

Mt
χ 2 (p)  k

k=1

σk2

,

(36)

normalized using adjustable weights σk that allow to rescale the relative contributions of each partial target

function to the total χ 2 . Naturally, it must hold that Mt > Mp in order to get meaningful results. Usual partial target functions are those that compare the time trace of the experimental FEB lineintegrated emissivities along each chord with those predicted by the numerical reconstruction. These functions are usually computed as integrals of the type ([t0 , tf ] denotes the temporal interval over which the reconstruction is done):

tf  = (χk2 )int(p)

 exp  dt  h(t  ) Lk (t  ) − LORCO (t  , p)  , k

t0

(37)

exp Lk (t  )

where is the experimental FEB time trace (t  , p)  is the time trace along the kth chord and LORCO k predicted by ORCO on the same chord using the current set of values for p.  LORCO is computed from k inserting the runaway density obtained from the transport model into Eq. (2) to get the local FEB emissivity:  = Cr ne (F, t)nr (F, t, p)Z  eff (t). IFEB (F, t, p)

(38)

Then, line-integration along the kth chord gives:

(t, p)  = IFEB (lk (F ), t, p)dl  k, LORCO (39) k kth-chord

where dlk is the length element along the relevant chord. Finally, the function h(t  ) in Eq. (37) is a temporal window that allows to focus the optimization effort within a prescribed temporal interval [td , tu ]. It is defined as: h(t) = H (t − td ) − H (t − tu ), t0
td < tu
tf ,

(40)

where H (t) is the usual Heaviside step function. Sometimes, it may also be important to reproduce the times at which the signals reach their maxima in intensity at the different chords. This can be done using as target functions the comparison of the times at which the maximum of the line-integrated FEB signals are reached: (χk2 )max (p)    exp   2 (t, p)  , = max Lk (t) − max LORCO k t ∈[t0 ,tf ]

t ∈[t0 ,tf ]

(41)

R. Sánchez et al. / Computer Physics Communications 156 (2003) 95–107

which can be useful to follow the propagation of these maxima as a runaway population generated at the central locations is transported outwards. Coming now to the free parameter vector, we have chosen as its first component the unknown constant appearing in Eq. (2): p1 = Cr . The rest of the components are provided by the representations of both D and VcF , that have been expressed in terms of expansions in orthogonal polynomials [31]: D(F, p)  =

lD


pl+1 Hl (F )

(42)

l=1

and VcF (F, p)  =

lV


pl+lD +1 H2l−1 (F ).

(43)

l=1

Thus, the total number of free parameters in the optimization is Mp = lD +lV +1. Notice that only odd polynomials have been included in the expansion for VcF  to enforce that it vanishes at the magnetic axis (see Section 4). Several choices for the polynomials Hl are available inside ORCO including the Legendre, Gegenbauer and Tchebyschev families.

103

7. Test case: JET limiter discharge 29586 As an example of application we have estimated D(F ) during the current-up phase of a low-density limiter JET discharge (#29586). The main parameters of this shot are B = 2.5 T, a = 1 m, b = 1.3 m, toroidal current Ip = 1–1.6 MA, line-integrated central density ne = (0.5–1.5) · 1019 m−3 , Te = 2–3 KeV and Zeff = 1.5–3.5. (Detailed temporal traces of the plasma current, loop voltage and electron density, as well as many other details can be found in Ref. [11].) The evolution of the electric field has been reconstructed by means of the JETTO code [24]. The main reasons for choosing this case for analysis are two. First, because the geometry of limiter discharges resembles more closely the axis-shifted elliptical geometry assumed by ORCO. Second, because one of the averaged runaway diffusion coefficient estimations mentioned in the introduction (see Section 1) was obtained with this JET discharge: Dr  0.2 m2 /s [11]. Therefore, the results obtained with ORCO can be compared with an independent assessment. The evolution of the experimental line-integrated FEB emissivity signals is shown in Fig. 3 for repre-

Fig. 3. Temporal evolution of the FEB emissivities (in a.u.) detected in vertical (left) and horizontal (right) chords.

104

R. Sánchez et al. / Computer Physics Communications 156 (2003) 95–107

sentative vertical and horizontal chords (the physical location of the chords can be seen in the sketch of the diagnostic shown in Fig. 1). The signal corresponds to the [133, 200] KeV energy window, since that is the one with best statistics. Clearly, the maximum of the signals is first reached at the central chords (see horizontal chord #5 and vertical chord #15). Correlating the time of appearance of these first maxima with the current and density time evolution during the discharge, it is realized that the runaway generation takes place in the central region during the current rampup phase, when the parallel electric field is strong and the density low. Their energy at that moment is estimated to be (1–2) MeV [11], and seem to reach quickly their radiation-limit energy at approximately 2.5 MeV [32]. It is also observed that the FEB signals from the outer chords reach their maxima at later times, suggesting the these electrons are transported outwards with time. In the analysis of Ref. [11] it was found that diffusive transport might be sufficient to explain the time evolution of the FEB signal at the central chords (in particular, only chords 15–18 were used, and a testand-error method was used to determine Dr ), whilst it failed to account for the behavior at the ones whose dominant contribution was supplied by the edge. To check this conclusion, we first use ORCO with only the central chords (vertical: #13–17; horizontal: #3– 7), imposing that VFc  = 0. It will however prove useful later, as a way to check the robustness of the results of the optimization. Coming now to the results of this first optimization, it is shown in Fig. 4 that the shape of D obtained depends at first heavily on the number of polynomials included in the expansion, k. Indeed, the total target function χ 2 decreases with increasing k until k  4 (see inset of the same figure). Then, the shape (and the value of χ 2 ) stays roughly the same until k > 8, beyond which convergence becomes more dubious due to the similarity between the number of dependent and independent variables. The resulting diffusion coefficient is not constant, but increases towards the edge, being able of simulating the FEB traces pretty well (see Fig. 5). As a test to confirm that the diffusive transport channel provides indeed a plausible dominant mechanism, we have repeated the calculations allowing a non-zero drift velocity profile. But, the best fit to the FEB signals is again obtained by [almost] the same

Fig. 4. Profiles for D that yield a minimum of the total target function for different number of polynomials (k) using only the central chords in discharge 29586. Inset: value of the total target function as a function of k.

shape for D and a negligible drift velocity. This insensitivity to the existence of a possible drift might be interpreted as a confirmation of the dominance of radial diffusion for the runaway core transport. To continue with the analysis, all vertical and horizontal chords are now included in the optimization. Again, we will set first VcF  = 0. The shape of D providing with a minimum for the total target function is shown in Fig. 6. Clearly, diffusivity remains similar to that obtained with the central chords (as it should be if it is to reproduce the traces at those chords). However, the value of χ 2 deteriorates notably since the total penalty function increases by a factor close to 8–10, even after normalizing to account for the extra number of chords. The reason must be looked for in the high value of the partial targets associated to the FEB signals of the chords #1, #9, #11 and #19 (see Fig. 7), that the algorithm reconstruct much more poorly: these are precisely the ones looking at the plasma edge! If we now repeat the same robustness test carried out before by allowing VcF  = 0, a non-negligible [outward] drift velocity profile is obtained, particularly at the plasma edge (see Fig. 8). The total partial target χ 2 is reduced around a 30% with the change, while D

R. Sánchez et al. / Computer Physics Communications 156 (2003) 95–107

105

Fig. 5. Comparison of experimental and simulated FEB emissivities at central chords for the diffusive coefficient shown in Fig. 4 with k = 6.

Fig. 6. Profiles for D that yield a minimum of the total target function as a function of k using all FEB chords. Inset: value of the total target function.

Fig. 7. Contributions to the total target function of the different FEB channels for the optimization done with k = 6 and all FEB chords.

changes appreciably. These results do not necessarily imply the physical meaningfulness for such a drift, but they suggest instead that diffusive transport may not be sufficient to explain the contribution to the FEB data

from edge runaways or that those edge signals are too much distorted by noise. Finally, we will compare these results with those reported in Ref. [11]. It is first necessary to convert our

106

R. Sánchez et al. / Computer Physics Communications 156 (2003) 95–107

Fig. 8. Test of robustness for the results obtained with all chords and k = 6. The profile of the drift velocity that yields the best results is plotted with a black thick line.

diffusive coefficient to real-space. This can be easily done by rewriting the argument F as a function of the cylindrical coordinates:   D = D F (R, Z) .

(44)

In Fig. 9, we compare the constant diffusive coefficient reported in Ref. [11] with the diffusive coefficient obtained previously, using k = 6, for the central chords. Three curves for D have been included, that respectively correspond to the values along the Z-axis and the R-axis. Notice that there are two different functions for the latter, since the plasma Shafranov shift prevents D to be symmetric with respect to the magnetic axis. Two conclusions can be extracted from the figure: (1) that diffusivity seems to be much larger (almost seven times) closer to the plasma edge (but it must be said that accuracy is worst there) and (2) that the values of D found at the low-field side of the core (that correspond to the region explored by the chords #15–18 that were used in that reference), ∼ 0.3 m2 /s are in good agreement with the previous estimation ∼ 0.2 m2 /s. Especially, when considering the experimental uncertainties, the cylindrical geometry and the absence of any optimization of the previous analysis.

Fig. 9. Profiles for D along the R-axis and Z-axis corresponding to the results obtained with k = 6 and using only the central chords. Dr obtained in Ref. [11] is also included (dashed line) as a reference.

8. Conclusions We have developed a code that can be used to estimate the spatial structure of the runaway radial diffusive coefficient in tokamaks. The code relies on the optimization of an internal transport model to best reconstruct some set of experimental measurements related to runaway electrons. It represents an improvement with respect to other codes due to the incorporation to the optimization effort of details of the plasma evolution through the discharge. In its current implementation, only discharges with an elliptical plasma shape are amenable of analysis. However, the main algorithm does not assume any type of geometry. Therefore, it could also be applied to discharges with different shapes (for instance, to D-shaped plasmas as those characteristic of divertor operation) if a different function F is chosen, or even if a fully numerical description of the equilibrium (such as that provided by EFIT [18] or other similar codes) is implemented. This upgrade would not present major difficulties thanks to the highly moduler structure of ORCO. Similar comments can be made about the diagnostic system used to supply to ORCO with experimental data on which to base the optimization. In the present version, a FEB

R. Sánchez et al. / Computer Physics Communications 156 (2003) 95–107

diagnostic system that measures line-integrated emissivities has been used. This will prove useful for the analysis of the data obtained with the FEB diagnostics currently being upgraded at FTU. But it would also be easy to modify the diagnostic module in ORCO to accommodate other runaway diagnostics, such as the infrared cameras existent in the Toroidal EXperiment for Technically Oriented Research (TEXTOR) [33]. Finally, we would like to make some comments about the possibility of obtaining some information of the velocity dependence of the runaway diffusive coefficient within the framework presented in this paper. Clearly, the possibility of introducing some v dependence on the transport coefficients, and of evolving a transport equation for f (r , v ) (the v⊥ averaged runaway distribution function instead of nr (r )) appears to be feasible. A similar optimization might then yield information, not only of the spatial dependences of D, but on its velocity dependences as well. This would allow not only to estimate turbulence levels, but also to test experimentally the reliability of Eq. (1). However, this approach is presently hopeless unless the experimental diagnostic system is able to collect separate information from runaways with different energies. The only information of this kind at our disposal is provided by the energy windows of the FEB diagnostic. But, in the case of the JET discharge we analyzed in the present paper, these windows cover only energies between 133–400 KeV, while the energy of the runaways seems to reach the radiation-limit at approximately 2.5 MeV as they diffuse out of the tokamak. An improved analysis could be available at JET in the near future, as the JET FEB cameras have been upgraded since the measurements described in these paper were performed: at present [34], the FEB energy range extends up to 6 MeV (with still four energy windows available) and some study of the velocity dependence of the runaway transport coefficient should be possible. As a final comment it is fair to say that, even with this lack of energy resolution, the results herein reported are by no means worthless. In many experimental cases similar to the one analyzed in this paper, independent diagnostics strongly suggest that the runaway population appears to be forming a quasi-monoenergetic beam during most of the analysis time [32]. In this sense, the value of D obtained by ORCO pertains to runaways with a well defined energy.

107

Acknowledgements Valuable discussions with M. Varela and advice on the Levenberg-Marquardt algorithm from S.P. Hirshman are acknowledged. Research supported by Spanish DGES Project No. FTN2000-0965.

References [1] H. Dreicer, Phys. Rev. 117 (1960) 329. [2] J.R. Martín-Solís, J.D. Alvarez, R. Sánchez, B. Esposito, Phys. Plasmas 5 (1998) 2370. [3] R. Jayakumar, et al., Phys. Lett. A 172 (1993) 447. [4] J.W. Connor, Plasma Phys. Contr. Fus. 35 (1993) B293. [5] A. Rechester, M. Rosenbluth, Phys. Rev. Lett. 40 (1978) 38. [6] H.E. Mynick, J.D. Strachan, Phys. Fluids 24 (1981) 695. [7] F.L. Hinton, R.D. Hazeltine, Rev. Mod. Phys. 48 (1976) 239. [8] C.C. Hegna, J.D. Callen, Phys. Fluids B 5 (1993) 1804. [9] R.W. Harwey, et al., Phys. Plasmas 7 (2000) 4590. [10] P. Helander, et al., Phys. Plasmas 7 (2000) 4106. [11] B. Esposito, et al., Plasma Phys. Contr. Fus. 38 (1996) 2035. [12] R. Jaspers, et al., Nucl. Fusion 33 (1993) 1775. [13] R.J. Bickerton, Plasma Phys. Control. Fus. 39 (1997) 339. [14] P.H. Rebut, R.J. Bickerton, B.E. Keen, Nucl. Fusion 25 (1985) 1011. [15] B. Esposito, et al., Phys. Plasmas 10 (2002) 2350. [16] W.D. D’Haeseleer, et al., Flux Coordinates and Magnetic Field Structure, Springer-Verlag, New York, 1991. [17] J. Freidberg, Ideal Magnetohydrodynamics, Plenum Press, New York, 1987. [18] L.L. Lao, et al., Nucl. Fusion 25 (1985) 1611. [19] P. Froissard, Ph. Thesis, Universite de Provence, AixMarseille, 1992. [20] R.M. Kulsrud, et al., Phys. Rev. Lett. 31 (1973) 390. [21] J.W. Connor, R.J. Hastie, Nucl. Fusion 15 (1975) 415. [22] S.V. Putvinski, P. Barabaschi, N. Fujisawa, N. Putvinskaya, M.N. Rosenbluth, J. Wesley, Plasma Phys. Control. Fusion 39 (1997) B157. [23] M.N. Rosenbluth, S.V. Putvinski, Nucl. Fusion 37 (1997) 1355. [24] G. Cenacchi, A. Taroni, Rep. JET-IR(88) 03, JET Joint Undertaking, Abingdon, 1988. [25] A. Airoldi, G. Cenacchi, Nucl. Fusion 37 (1997) 1117. [26] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes, Cambridge Univ. Press, Cambridge, 1986. [27] F. Gill, Practical Optimization, Prentice-Hall, New York, 1989. [28] S.P. Hirshman, et al., Phys. Plasmas 6 (1999) 1858. [29] R. Sanchez, S.P. Hirshman, A.S. Ware, L.A. Berry, D.A. Spong, Plasma Phys. Contr. Fusion 42 (2000) 641. [30] A.H. Reiman, et al., Phys. Plasmas 8 (2001) 2083. [31] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions, Natl. Bureau of Standards, Washington, DC, 1968. [32] J.R. Martín-Solís, B. Esposito, R. Sánchez, J.D. Alvarez, Phys. Plasmas 6 (1999) 238. [33] R. Jaspers, Ph. Thesis, Technical University of Eindhoven, 1995. [34] V.G. Kiptily, Private communication, 2003.