Gyrokinetic full-torus simulations of ohmic tokamak plasmas in circular limiter configuration

Computer Physics Communications 203 (2016) 128–137

Contents lists available at ScienceDirect

Computer Physics Communications journal homepage: www.elsevier.com/locate/cpc

Gyrokinetic full-torus simulations of ohmic tokamak plasmas in circular limiter configuration T. Korpilo a,∗ , A.D. Gurchenko b , E.Z. Gusakov b , J.A. Heikkinen e , S.J. Janhunen c , T.P. Kiviniemi a , S. Leerink a , P. Niskala a , A.A. Perevalov d a

Department of Applied Physics, Aalto University, Espoo, Finland

b

Ioffe Physical-Technical Institute of the Russian Academy of Sciences, St. Petersburg, Russia

c

Princeton Plasma Physics Laboratory, Princeton University, NJ, USA

d

St. Petersburg State Polytechnic University, St. Petersburg, Russia

e

VTT, Espoo, Finland

article

info

Article history: Received 6 March 2015 Received in revised form 27 January 2016 Accepted 22 February 2016 Available online 7 March 2016 Keywords: Core–edge turbulence Gyrokinetic particle simulation Fusion plasma physics

abstract The gyrokinetic full 5D particle distribution code ELMFIRE has been extended to simulate circular tokamak plasmas from the magnetic axis to the limiter scrape-off-layer. The predictive power of the code in the fulltorus configuration is tested via its ability to reproduce experimental steady-state profiles in FT-2 ohmic L-mode plasmas. The results show that the experimental profile solution is not reproduced numerically due to the difficulty of obtaining global power balance. This is verified by cross-comparison of ELMFIRE code versions, which shows also the impact of boundary conditions and grid resolution on turbulent transport. © 2016 Elsevier B.V. All rights reserved.

1. Introduction Direct numerical simulations based on gyrokinetic theory have become an important tool for studying turbulent transport physics in tokamak plasmas. Computation power of modern computers together with state-of-the-art numerical techniques enables addressing the increased complexity of such simulations, characterized by transitions from δ f to full f method, from adiabatic to kinetic electrons, from local to global model, from analytic magnetic background to full magnetic geometry and from electrostatics to electromagnetism. Despite these advances, gyrokinetic codes are still mostly restricted to the core plasma region with closed magnetic flux surfaces. To capture the complex picture of plasma dynamics in a tokamak, coupled core–edge simulations are however necessary because the edge plasma and the material boundary determine the conditions in the core plasma, while cross-field transport from the core determines the conditions in the edge plasma. Moreover, as the edge plasma with steep radial gradients and large-amplitude fluctuations is not in thermodynamic equilibrium, full electron and ion distribution function, full f , kinetic ap-

∗

Corresponding author. E-mail address: [email protected] (T. Korpilo).

http://dx.doi.org/10.1016/j.cpc.2016.02.021 0010-4655/© 2016 Elsevier B.V. All rights reserved.

proach is required in such simulations. As an effort towards better understanding of tokamak plasma transport, kinetic and fluid codes are being developed for core–edge plasma turbulence [1,2]. To include the core–edge interaction, the electrostatic full f gyrokinetic particle-in-cell code ELMFIRE [3] has been extended to simulate tokamak plasmas from the magnetic axis to the material boundary in an axisymmetric toroidal configuration. The material boundary consists of a radial wall and discrete poloidal limiter plates, the latter of which define the scrape-off-layer and the last closed flux surface. As a physical simplification of the simulation model, the thin net-charge layer at the plasma–material interface is omitted from the self-consistent solution. This is justified because the Laplace operator, that would enable the charge separation, is omitted from the Poisson equation, as commonly done in gyrokinetic realizations. The sheath potential then arises from the quasineutrality condition, and locally high grid resolution is not required to resolve the thin Debye sheath. Furthermore, the important edge plasma elements of neutral transport and atomic processes are not modeled. Several ELMFIRE predictions on core turbulent transport phenomena have been recently validated against measurements from the FT-2 experiment [4]. Despite this, the simulations have not reproduced the experimental level of transport required for steady-state density and temperature profiles. However, as conducted in a toroidal annulus with uniform grid spacing, the

T. Korpilo et al. / Computer Physics Communications 203 (2016) 128–137

simulations have been under-resolved in grid resolution while sources and sinks have been modeled inaccurately at the radial boundaries. The impact of these inaccuracies on the transport results has been an open issue but can now be addressed with the full-torus version of the code. The present full-torus simulations enable resolving the plasma turbulence up to the local ion gyroradius scale across the poloidal plane and include more accurately the sources and sinks, in particular the sink action of the limiter plates. In the present paper, the numerical methods developed to include the material boundary and magnetic axis region in ELMFIRE are presented. The predictive capability of the code in the resulting full-torus configuration is tested by simulating ohmically heated FT-2 limiter plasma with the aim of obtaining a numerical steady-state solution close to the experimental profile equilibrium when both the heat and particle sinks and sources are taken into account. The impact of boundary conditions and grid resolution on particle and heat transport is shown in the process. The paper is organized as follows. In Section 2, the physical model and numerical methods used in the ELMFIRE code are described. In Section 3, the implementation of magnetic axis region, radial wall and limiter plates into the code is presented. In Section 4, the numerical reproduction of experimentally measured steady-state profiles is attempted. In Section 5, a cross-comparison of ELMFIRE code versions is performed to verify the obtained results. Summary and conclusions are given in Section 6. 2. ELMFIRE simulation model The gyrokinetic equations solved in the ELMFIRE code are based on the second-order gyrokinetic formalism in which the effect of the ⟨E⟩ × B drift velocity is eliminated from the definition of the gyroradius [5]. This introduces the polarization drift in the equations of motion and eliminates the explicit appearance of the polarization density operator from the gyrokinetic Poisson equation. The gyrokinetic equations of motion for the gyrocenter position R and parallel velocity U are given in the formalism as

˙ = R U˙ =

ˆ + Eˆ × b BU B⋆ ˆ · Eˆ q B m

B⋆

ε0 ∇ · E (r, t ) =

  q

f δ r − x(2)

 



d3 RdUdJ ,

(3)

p.s.

which is the quasineutrality condition — the divergence of the electric field is neglected in numerical realizations. Here, the summation is over the particle species, ε0 is the vacuum permittivity, δ is the delta function, and f (R, U , J , t ) denotes the primary form of the distribution function that is resolved in particle simulations. In the present paper, the long wavelength limit of perturbations is assumed as it simplifies the gyrokinetic equations by making the second-order correction ∆ρ in x(2) to vanish. The ELMFIRE code simulates the full particle distribution function of drift-kinetic electrons and gyrokinetic ions in the presence of a 3D electrostatic potential in an axisymmetric toroidal magnetic equilibrium with circular and concentric flux surfaces. The straight field line coordinate system in such a magnetic equilibrium takes the form r = (r , θ , φ), where the coordinates denote radial coordinate, poloidal-like angle and toroidal angle, respectively. The equations of motion are written in terms of the straight field line coordinates, thus resembling the Hamiltonian form, and time integrated by a hybrid implicit–explicit particle mover. The implicit particle mover in the integration scheme is used to neutralize the charge separation between ions and electrons caused by the explicit particle push. From the terms present in the equations of motion, the ion polarization drift velocity and the electron acceleration by the parallel electric field (and the ion acceleration for symmetry) are integrated implicitly in the code; the former is integrated by the implicit Euler method and the latter by the implicit Velocity Verlet method [6]. The rest of the terms in the equations of motion are integrated explicitly by the fourth-order Runge–Kutta method, RK4. The algorithm for the ion particle mover is then

˙ − vpol (zn ) Rn+1 = Rn + RK4 R 



1

+ a∥ (zn+1 ) ∆t 2 b (Rn+1 ) + vpol (zn+1 ) ∆t ,

,

Un+1 = Un + RK4 U˙ − a∥ (zn ) +



,

ˆ · b and where the fields Eˆ and Bˆ are defined as where B⋆ = B J 2π m 

∇B −

m d

E x(1)

 

q dt



B

  (1)  E x ×b ˆB = B + m ∇ × Ub + . q

ˆ The drift-kinetic, limit and by omitting the last term of Eˆ and B. gyrokinetic Poisson equation is given in the adopted formalism as

(4)

2

(1)

   Eˆ = E x(2) −

129



×b

 , (2)

B

Here, J = π mV⊥2 /Ω is the adiabatic invariant (J˙ = 0), ⟨. . .⟩ operator denotes the averaging with respect to the gyroangle, b is the unit vector parallel to the magnetic field B = Bb, E is the electric field and Ω = qB/m is the gyrofrequency, while V⊥ , q and m denote perpendicular velocity, charge and mass, respectively. The magnetic field related variables are defined at the gyrocenter position R while, according to the gyrokinetic ordering, the electric field is averaged over the positions x(1) = R + b × V⊥ /Ω or x(2) = x(1) + ∆ρ. The term ∆ρ consists of short wavelength trembling that remains in the gyroradius after the elimination of the ⟨E⟩ × B drift effect. The magnetic drift terms arising from an inhomogeneous ambient magnetic field are included in the above equations of motion as they appear in the standard gyrokinetic theory. The standard guiding center equations of motion can be reduced from Eqs. (1) and (2) by taking the zero-gyroradius, or



1 2

a∥ (zn+1 ) + a∥ (zn ) ∆t ,



where vpol = (m/qB⋆ ) b × d (⟨E⟩ × b/B) /dt is the polarization drift velocity, a∥ = (q/mB⋆ ) B · ⟨E⟩ the acceleration by the parallel electric field, ∆t the time step size, z represents the phase-space coordinates (R, U , J ), and the subscript n indicates the time level. The convective time derivative of the ⟨E⟩ × B drift term is discretized by the backward difference scheme in the polarization drift. The algorithm for the electron particle mover is identical to the ions but ˙ and because the standard guiding center equations of motion for R U˙ are solved for electrons the polarization drift is not present in the algorithm. The quasineutrality condition and the related electric potential are solved in a structured toroidal grid set in the straight field line coordinates. The grid supports a nonuniform spacing between the radial grid points and a radially varying poloidal grid size, thus allowing grid resolution of the order of local ion gyroradius across the poloidal cross section of the simulated plasma. The interpolation between the grid and  a particle is carried out by a trilinear shape function S rp − rk where rp and rk denote the particle and grid point position, respectively.   Here, the electron position is determined by the point δ rp − R , while four points

on the gyroring δ rp − x(2) are used to describe an ion. The linear interpolations in the shape function are performed along the straight field line coordinates (r , θ) and along the magnetic field, instead of the toroidal direction, in order to make use of the long parallel wavelength characteristic of perturbations in the core





130

T. Korpilo et al. / Computer Physics Communications 203 (2016) 128–137

plasma. The magnetic shear is recovered correctly in the process when the interpolation in the radial direction is performed first. The covariant components of the electric field are calculated by a momentum conserving scheme in the explicit ⟨E⟩ × B drift evaluation and by an energy conserving scheme in the implicit parallel acceleration and ion polarization     Both  drift evaluation. schemes take here a form of E rp = E (rk ) S rp − rk where the electric field E (rk ) is calculated from the potential grid values ϕ (rk ) by the central difference method. The difference between the momentum and energy conserving schemes is that, in the former, the electric field E (rk ) is defined at the grid point whereas, in the latter, at the center between two grid points. For numerical stability, the formula for the radial electric field in the momentum conserving scheme is modified so that the condition (∇ × E)φ = 0 is satisfied [7]. The electric potential is solved in the code from the quasineutrality condition ρe (rk ) + ρi (rk ) = 0 by means of the implicit particle push, cf. Eq. (4). The electron charge for example, is   density, calculated on the grid as ρe (rk ) = qp wp S rp − rk /Vk where wp and Vk denote the particle weight and grid cell volume, respectively. Writing the electron quasineutral position as rp = r′p + drp , where r′p is the electron position after the Runge–Kutta step and drp the small displacement due to the acceleration by the parallel electric field, the Taylor series expansion of the shape function up to the first order term gives the electron charge density, and in similar manner the ion charge density, as

ρe (rk ) = ρ¯ e (rk ) +



drp z′p ·

 

p

ρi (rk ) = ρ¯ i (rk ) +





qp wp ∂ S (r − rk ) Vk

∂ζ

 ∇ζ r=r′p

,

dsp z′p

 

p

  qp wj  ∂ S (r − rk ) ∂ S ( r − rk ) · ∇r + ∇θ . Vk ∂r ∂θ r=r′ j

(5)

p,j

Here, the coordinate ζ denotes the angle along the magnetic field line, j indexes the four points on the ion gyroring, drp =   1 qE∥,n+1 /m ∆t 2 b and dsp = vpol ∆t ≈ (⟨E⊥ ⟩n+1 − ⟨E⊥ ⟩n ) /Ω B. 2 The charge densities ρ¯ i and ρ¯ e are evaluated at particle positions r′p,j and r′p , respectively. The small effect of the ion acceleration by the parallel electric field is excluded from the displacement dsp , and the electric field perpendicular to the magnetic field, E⊥ , is approximated as a radial and poloidal component in the polarization drift. The quasineutrality condition together with Eq. (5) gives a linear system for the electric potential ϕn+1 , from which the potential is solved by conventional means. After solving the electric potential, the implicit terms in Eq. (4) are evaluated. Collisions between charged particles are modeled in the driftkinetic limit by the binary collision model introduced in Refs. [8,9]. The total toroidal ohmic current is induced in the code by a time varying, spatially homogeneous loop voltage [10], but instead of using the self-consistent toroidal current, the poloidal magnetic field is calculated from an input toroidal current profile. The conservation of total energy and total toroidal angular momentum in ELMFIRE simulations is reviewed in Refs. [6,7]. 3. Numerical realization of material boundary and magnetic axis region The ELMFIRE simulation domain covers the entire tokamak plasma volume from the magnetic axis to the material boundary of a radial wall and poloidal limiter plates. Numerically, the merit of incorporating the plasma center into the gyrokinetic simulation lies in the absence of the inner radial simulation boundary which not only removes the stability issues of that boundary but also

makes the remaining simulation boundary more stable. However, the extension of the simulation domain to the magnetic axis introduces a coordinate singularity, while the extension to the material surface includes the formulation of boundary conditions that are suitable for a gyrokinetic particle approach. The magnetic axis region is included in the ELMFIRE simulation model by using the straight field line coordinates. To avoid resulting grid singularity, the simulation grid is taken to the plasma center by decreasing the poloidal grid size radially so that there is only one poloidal grid point at the magnetic axis (r = 0). The interpolation between the grid and a particle is performed by the trilinear shape function S also at the axis, although the grid structure is triangular there. In addition to the grid, the singularity in the straight field line coordinates affects the time integration of the particle motion because the poloidal equation of motion and the poloidal electric field are proportional to r −1 . To overcome the impact of the coordinate singularity on the particle motion, the time step is shortened in the explicit Runge–Kutta integration for particles affected by the singularity, and particles that get closer than one micrometer to the magnetic axis are moved over the axis. Considering the implicit particle motion, the acceleration by the parallel electric field and the concomitant displacement drp are not singular, while the poloidal component of the polarization shift dsp is. The singularity in the polarization shift does not, however, affect the gyrokinetic Poisson equation because the inverse of the radial coordinate appears in the equation together with the poloidal derivative of the shape function in the form of r −1 ∂ S /∂θ . The cancellation of the r −1 effect near the magnetic axis is then explained by writing out the shape function ∂ S /∂θ = (1 − r /∆r ) (∂ S /∂θ)0 + (r /∆r ) (∂ S /∂θ)1 in terms of the radial interpolation for the interval 0 ≤ r ≤ r1 . Here, ∆r is the length of the interval, the subscripts 0 and 1 refer to the radial grid points at the magnetic axis and the next, respectively, and (∂ S /∂θ)0 = 0. The singularity of the poloidal electric field at the magnetic axis is avoided for the same reason. Another complication arising from the coordinate singularity concerns the basis vectors of the straight field line coordinates because their direction and magnitude may change significantly between the points on the gyroring and the gyrocenter near the magnetic axis. Therefore, in the evaluation of the gyroaverage of the electric field, the radial and poloidal basis vectors at the gyroring points are first transformed to the basis vectors of the standard cylindrical coordinates (R, Z , φ), and after the averaging process the cylindrical basis vectors are transformed back to the radial and poloidal basis vectors at the gyrocenter position. A similar procedure is performed in the gyrokinetic Poisson equation when taking the dot product of the gyrocenter polarization shift and the radial and poloidal basis vectors at the gyroring points. A wall surface at fixed radius parallel to the magnetic field is taken as the radial simulation boundary in the code. The simulation grid is restricted radially to the wall in such a way that the first radial grid point omitted in the calculation is at the wall surface. To avoid large charge imbalance near the wall boundary and related numerical instabilities, the gyrocenter orbits of ions and the guiding center orbits of electrons are followed and kept in the simulation until they stop contributing charge to the simulation grid. Ions are thus lost to the wall when their entire gyroring is inside the wall, while electrons are lost to the wall when their guiding center crosses the wall surface. As the potential boundary condition, zero electric potential is applied at the wall surface and inside the wall. However, because the polarization drift is proportional to the time derivative of the electric field, the boundary condition of zero electric potential has only a weak impact on the drift and, therefore, on the gyrokinetic Poisson equation. To enhance the effect of the potential boundary condition in the Poisson equation, the electric field term ⟨E⊥ ⟩n in

T. Korpilo et al. / Computer Physics Communications 203 (2016) 128–137

the polarization shift dsp is set artificially to zero for ions whose gyrocenter is inside the wall. A toroidal limiter is modeled in the code effectively as a toroidal chain of identical poloidal limiters as illustrated in Fig. 1. The advantage of poloidal limiters over a continuous toroidal limiter is the compatibility with the applied shape function, which interpolates along the magnetic field in one direction, making poloidal limiters numerically more feasible. Here, the poloidal limiters are taken as toroidally infinitely thin plates that are limited radially between the last closed flux surface and the wall and poloidally between the angles θ1 and θ2 . The poloidal size of the plates is determined so that the magnetic field lines inside the scrape-off-layer region cannot pass between two consecutive plates without hitting one of them. The limiter plates are included in the simulation by placing them at the positions of the toroidal grid points. On these poloidal planes, the plate boundaries lie on the grid lines following the poloidal grid points at the last closed flux surface and the radial grid points at the angles θ1 and θ2 . Electrons (ions) are considered lost to the limiter plate when the guiding center (any of the four points on the gyroring) hits the plate. Consequently, because the grid points on the limiter plates and their boundaries are omitted in the calculation, the electron charge flowing to the limiter plates vanishes gradually from the simulation grid, while sudden ion charge losses occur at the grid points near the limiter plates. As the potential boundary condition, the limiter plates and their boundaries are set to zero electric potential. The effect of the potential boundary condition on the gyrokinetic Poisson equation is enhanced again by the artificial manipulation of the electric field in the polarization drift. In other words, if a point (out of four) on the ion gyroring lies closer than a certain poloidal distance from the poloidal boundaries of a plate when its position is projected on the poloidal plane of the plate, the electric field from that poloidal plane is zeroed for the point in question in the evaluation of the gyroaverage ⟨E⊥ ⟩n . The certain poloidal distance at each radius is determined adaptively so that the potential grid values next to the poloidal boundaries of the plates remain at the average grid potential in front of the plate surfaces, and it is equal to one poloidal grid cell width at most. For numerical stability, the ⟨E⟩ × B drift velocity is suppressed artificially near the limiter plates by zeroing the electric potential at the poloidal grid points next to the plate poloidal boundaries when evaluating the drift for both ions and electrons. 4. Comparison of ELMFIRE results with experimental data The predictive capability of full-torus ELMFIRE simulations in the limiter configuration is investigated here by testing the code’s ability to reproduce measured steady-state profiles. In this test, a simulation is started from experimental equilibrium solution for simplicity and then evolved self-consistently in time in the presence of heat and particle sources and sinks. For the experimental profiles to be maintained in the simulation, the main plasma processes need to be modeled accurately and the effect of the sources and sinks has to balance out. Naturally, the parameters of the experimental case also have to be suitable for ELMFIRE in order to have optimal conditions for the simulation–experiment comparison. The research FT-2 limiter tokamak [11] – a large-aspect-ratio, low-beta machine with a circular plasma cross section – is chosen here as the testbed for studying the ELMFIRE code. An advantage of the tokamak is large ρ∗ plasmas that are computationally inexpensive for a gyrokinetic full f particle-in-cell code to simulate with grid resolution of the order of ion gyroradius across the entire poloidal cross section. Another advantage is that the fluctuation levels are high in FT-2 which relaxes the number of particles required to overcome noise in the simulation. To have then the

131

Fig. 1. The material boundary in ELMFIRE simulations: toroidal limiter is approximated by a toroidal chain of identical poloidal limiters (in red) protruding from the wall (in blue). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

optimum experimental plasma equilibrium and conditions for the simulation, a tailor-made ohmic L-mode FT-2 discharge has been executed. In this discharge, the magnetic geometry is characterized by circular flux surfaces centered close to the magnetic axis, making the approximation of concentric circular flux surfaces in the code valid. The electrostatic approximation, in turn, is justified because of the low-beta plasma and because the observed MHD activity is low in the discharge. Considering the FT-2 limiter structure, two poloidal ring limiters are located toroidally at the opposite side of the tokamak vessel. The toroidal limiter model in Section 3 is, however, used in the simulation because the real poloidal limiter configuration is found to suffer from numerical problems. To match the plasma-limiter contact point, the limiter plates are placed at the high-field-side midplane in the simulation as the radially outer magnetic flux surfaces are shifted slightly inwards in the discharge. The tailor-made FT-2 hydrogen discharge is characterized by low current, high density and high collisionality, ν ∗ = 5 − 30. The main parameters of the discharge are major radius 55 cm, limiter radius 7.9 cm, on-axis magnetic field 2.0 T and toroidal plasma current 28.4 kA. The measured density and temperature profiles as well as the plasma current profile based on ASTRA [12] modeling are shown in Figs. 2–3 after mapping them to the concentric flux surfaces. The profile data points from 7.5 cm onwards are based on extrapolation and not on experimental information. The plasma impurity concentration parameter is observed to be Zeff = 2.2, and the impurity content is simulated as highly ionized oxygen O6+ ions as they are the dominant charge state in the discharge. The steady-state oxygen temperature is assumed to be equal to the hydrogen temperature, and the oxygen density is assumed to have the profile of the electron density. To adopt the experimental equilibrium in the simulation, the limiter radius is set to 7.6 cm which is approximately the average last-closed-flux-surface radius in the discharge. The wall boundary is placed at the radius of 8.1 cm, and the 5 mm thick scrape-off-layer is defined by 16 limiter plates with the poloidal size of ∆θ /2π = 0.03. Hence, the simulation domain covers the entire plasma volume from the magnetic axis to the material boundary of the radial wall and the limiter plates. The applied toroidal grid is composed of around 0.9 million points in 16 toroidal sections, which provides perpendicular resolution of proton gyroradius up to about 7 cm in radius after which the poloidal grid size is limited to 900 points and radial resolution to 0.4 mm. A time step ∆t = 15 ns is chosen for the simulation to satisfy the accuracy conditions of k∥ vte ∆t . 1 and k⊥ vE ×B ∆t . 1

132

T. Korpilo et al. / Computer Physics Communications 203 (2016) 128–137

Fig. 2. Measured electron and hydrogen density (left) and temperature profiles (right) along with initial and final ELMFIRE profiles.

Fig. 3. Steady-state ELMFIRE plasma current profile (left) and final ELMFIRE profiles of radial energy flux for electrons and hydrogen ions (right) along with ASTRA predictions.

Fig. 4. Measured radiation profile (left), and the applied probabilistic ionization rate distribution (right).

where k∥ denotes the largest wavenumber of the parallel electric field and vte the electron thermal velocity. The energy confinement time of FT-2 discharges is of the order of 1 ms. The total number of simulated electrons is 5.8 billion, translating into 30 000–500 (axis–edge) electrons per grid point. In ohmic FT-2 discharges, the energy balance is determined by ohmic heating, radiation, and energy loss to the limiters and chamber wall. From these sources and sinks, the ohmic current is

included in the simulation by an active feedback control of the loop voltage which ramps up and sustains the toroidal plasma current of 28.4 kA. The radiation losses, in turn, are modeled by cooling the electrons according to the experimentally measured radiation loss rate profile shown in Fig. 4. The energy loss to the limiter plates and the wall by particles is inherently incorporated into the simulation. The particles lost to the limiter plates and the wall are recycled back to the plasma as 2 eV neutral ion–electron pairs according to

T. Korpilo et al. / Computer Physics Communications 203 (2016) 128–137

133

Fig. 5. Time evolution of total energy (left) and total toroidal angular momentum plus cumulative momentum transport across the diagnostic domain boundary (right).

Fig. 6. Analysis of particle balance ∂t n + ∇ · Γr = S (left) and power balance ∂t (1.5nT ) + ∇ · Qr = P (right) for electrons and hydrogen ions inside the last closed flux surface at the end period of the simulation. Γr , Qr and u denote radial particle flux, radial energy flux and energy density, respectively. Particle source S comes from particle recycling, and the term P includes ohmic power, radiation power losses (electrons only) and power transfer through particle collisions.

the radial ionization rate profile in Fig. 4 that is deduced from a FT-2 neutral hydrogen density profile [13]. In the ELMFIRE simulation, the particles are initialized on their collisionless neoclassical orbits so that the initial density and temperature profiles match exactly the experimental ones. The particle distribution is then let to evolve self-consistently in the presence of the sources and sinks for 93 µs, a fraction of the energy confinement time, within which the turbulent and neoclassical transport levels are well saturated. The resulting plasma profiles are shown in Figs. 2–3. In summary, the toroidal plasma current profile reaches a steady-state solution in the simulation that is in close agreement with the ASTRA current profile used to calculate the poloidal magnetic field in ELMFIRE. The loop voltage of 4.1 V is found to drive the 28.4 kA plasma current at the end of the simulation while the experimentally observed loop voltage is 3.7 V. As an indicator of the numerical quality of the simulation, the diagnosed toroidal angular momentum is conserved within 0.5% accuracy as seen in Fig. 5. The conservation is examined by recording the momentum inside the radius of 7.3 cm as well as the cumulative momentum transport through that radial boundary (the effect of radiation losses is small). The density profiles, in turn, are well maintained in the simulation albeit small profile relaxation at the plasma edge; the subtraction of the hydrogen and electron densities gives effectively the oxygen charge density profile. The quasi-steady density profiles are partly the result of reasonable compensation of the radial particle fluxes by the

recycling particle sources as shown by the particle balance analysis in Fig. 6. The simulation time is also very short compared to the particle confinement time, which is a few energy confinement times in FT-2 discharges. According to the diagnostics recording the plasma energy content inside the radius of 7.3 cm in Fig. 5, the main issue in the simulation is that the power balance is not achieved. In other words, a net heating power of 58 kW is observed inside the diagnostic domain, which results from ohmic heating power of 113 kW, numerical heating of 3 kW, radiation power losses of 16 kW, and transport losses of 42 kW across the domain boundary. The reason for the power imbalance is the underprediction of electron radial energy (convective plus conductive) transport as shown by the power balance analysis in Fig. 6. Transport analysis with the ASTRA code in Fig. 3 indicates also that much higher electron energy flux across the radius is required for the steadystate solution when the same power sources and sinks are applied in ASTRA and ELMFIRE. The power imbalance manifests itself as electron heating clearly visible in the electron temperature profile. The ion power balance analysis, in turn, shows an excessive source of radial energy transport at mid-radii and, like electrons, a transport sink at the outer core region. This mechanism causes the observed flattening of the hydrogen temperature profile; the oxygen temperature profile shows similar relaxation. The ASTRA modeling indicates that slower energy transport from r = 3 cm onwards would be required to maintain the measured hydrogen

134

T. Korpilo et al. / Computer Physics Communications 203 (2016) 128–137

Fig. 7. Neoclassical radial particle flux (left) and radial heat flux (right) profiles for electrons and hydrogen ions from ELMFIRE and NCLASS. The neoclassical ELMFIRE simulation is performed by filtering out non-zero poloidal and toroidal mode numbers in the term ρ¯ e + ρ¯ i in Eq. (5). The simulation domain is also limited to the last closed flux surface, and ohmic heating (loop voltage) and radiation losses are disabled.

Fig. 8. Electron density fluctuations on a poloidal cross section at t = 93 µs given by the present ELMFIRE version (left) and an older code version (right).

temperature profile without any additional energy sources and sinks. To check that the underlying neoclassical transport is modeled correctly in the simulation, the ELMFIRE particle and heat fluxes in the neoclassical limit are compared with the predictions of the NCLASS code [14] in Fig. 7. The comparison shows a decent agreement between the code results, considering that ion orbits are very wide relative to the macroscopic scale lengths in FT2 plasmas which affects the validity of the analytic neoclassical theory [15]. As for turbulent transport, the plasma fluctuations associated with anomalous energy and particle fluxes are shown in the left-hand panel of Fig. 8. The linear frequency and growth rate spectra of these fluctuations are also calculated in the collisionless limit and compared with the predictions of the gyrokinetic GS2 code [16] in Fig. 9. The linear analysis is performed here in two distinct parameter regimes at r = 3 cm and r = 6 cm, and the linearly unstable modes in these two regimes are unlikely overlapping as the radial envelope of the modes is bounded in r. The results show a fair agreement of linear frequencies and growth rates between ELMFIRE and GS2 at the radius r = 3 cm where the ion temperature gradient is mild. In contrast, at the radius r = 6 cm with steep ion temperature gradient, the frequency match is poor between the codes although the growth rate spectra are in fair agreement. The frequency match is, however, significantly improved by increasing the plasma current, which suggests that non-local effects modify the ELMFIRE frequencies.

Judging from the results, the simulation is maintaining neither the measured steady-state temperature profiles nor the power balance, and consequently the experimental profile solution is not reproduced numerically. An analogous result is obtained from a simulation of another tailor-made ohmic L-mode FT-2 discharge that has lower density and higher heat transport level than the discharge considered here. The reason for poor power balance or underprediction of electron energy transport across the radius in these simulations is not presently understood. One possible explanation is the missing sub-ion-scale turbulence [17]. By contrast, given the low beta and high density of the discharge, electromagnetic turbulence and runaway electrons are unlikely to explain the electron transport underprediction. It is equally unlikely that the turbulence drive would be underestimated at all radii in the simulations due to the experimental errors in the measured input parameters. The reason for the strong relaxation of the ion temperature profiles is presently also unclear but the heating of the outer core and plasma edge region could be attributed to the absence of neutrals and atomic and molecular processes, such as charge exchange reactions, or to the transport shortfall described in Ref. [18]. In addition, shorter time step, higher particle number and more accurate limiter model in particular can have an effect on the time evolution of the edge plasma parameters.

T. Korpilo et al. / Computer Physics Communications 203 (2016) 128–137

135

Fig. 9. Comparison of ELMFIRE and GS2 linear frequency and growth rate spectra of electric potential fluctuations at radii r = 3 cm (left) and r = 6 cm (right). The ELMFIRE results are obtained without filtering from the linear phase of the simulation that does not include ohmic heating, radiation losses, binary collisions and scrape-off-layer. ρH , vtH and R denote proton gyroradius, proton thermal velocity and major radius, respectively, and negative frequency is in the electron diamagnetic direction.

Fig. 10. Final electron and hydrogen density (left) and temperature profiles (right) in the full-torus and annulus simulation.

5. Cross-comparison of ELMFIRE simulations To verify the presented full-torus simulation results, the ohmic L-mode FT-2 experiment is re-examined by using an older version of the code. As the direct predecessor of the present version, the older code version implements the equations and algorithms described in Section 2, except for the axisymmetric toroidal grid. The older version uses instead a field-aligned (quasi-ballooning) grid structure that has a uniform radial grid spacing and a radiallyfixed poloidal grid size. Another significant difference between the code versions is that the older code version simulates a toroidal annulus between an inner and outer radial boundary, which means the plasma center and scrape-off-layer regions are not included in the simulation domain. The boundary treatments applied at the inner and outer simulation boundary, different from those discussed in Section 3, are described in Refs. [3,7]. As the potential boundary conditions, the outer radial boundary is set to zero potential while zero radial electric field is imposed at the inner radial boundary. In addition, particles crossing the outer radial boundary are recycled back to the plasma as cold neutral electron–ion pairs. The FT-2 case is rerun with the older code version by taking the discharge parameters and profiles described in the previous section as initial conditions for the simulation. The heat and particle sources and sinks are modeled as in the full-torus simulation, given that the outer boundary acts as the particle sink instead of the wall and the limiter plates. As for numerical

parameters, the simulation domain is limited between the inner radial boundary at 2 cm and the last closed flux surface at 7.6 cm. The simulation grid is composed of conventional 110 radial, 150 poloidal and 8 toroidal points, which gives fixed 0.5 mm radial resolution across the simulation domain. For comparison, the full-torus simulation provides 0.9–0.4 mm radial resolution and 145–900 poloidal grid points over the same region, enough to capture the ion-scale dynamics. The total number of simulated electrons is 0.35 billion, translating into 1000 particles per grid point at the outer core region, which is equal to the full-torus simulation. The time step of 30 ns is used in the simulation. The simulation with the older code version, referred to as the annulus simulation, is initialized as the full-torus simulation, and thus the initial plasma profiles match between the simulations. The final profiles are shown in Figs. 10–12, where the full-torus simulation profiles of density, temperature and plasma current are the same ones as in Figs. 2–3. On the basis of the results, the toroidal current profiles are similar in the simulations; the loop voltage in the annulus simulation is about 3.8 V (vs. 4.1 V) at the end of the simulation which is close to the measured value. The density profiles are also similar between the simulations although some particle accumulation is observed near the outer radial boundary in the annulus simulation. The toroidal angular momentum is conserved well in both simulations, within 0.6 % in the annulus simulation, while the diagnostics recording the total energy content in the annulus simulation shows an identical power imbalance issue to the full-torus simulation. In other words,

136

T. Korpilo et al. / Computer Physics Communications 203 (2016) 128–137

Fig. 11. Final plasma current profiles in the full-torus and annulus simulation (left), and time evolution of total energy in the annulus simulation (right).

Fig. 12. Final profiles of radial particle flux (left) and radial heat flux (right) for electrons and hydrogen ions in the full-torus and annulus simulation. Inside the last closed flux surface, the difference between electron and hydrogen radial particle flux is explained by oxygen flux not shown in the figure.

ohmic heating, radiation cooling and energy transport losses are not balanced out in either simulation, resulting in similar electron heating and evolution of the electron temperature profiles. The hydrogen temperature profiles show similar relaxation in both simulations. Judging from the results, the power imbalance and associated electron heating as well as the relaxation of the ion temperature profiles are reproduced with the older code version, which verifies the earlier results of the simulation–experiment comparison. In addition to this, the comparison of the two code versions reveals the impact of the boundary conditions and grid resolution on the simulation results. First, Fig. 12 shows that the radial particle fluxes are strongly reduced towards the outer simulation boundary in the annulus simulation, which explains the observed particle accumulation. In the full-torus simulation, by contrast, the particle fluxes are much higher at the outer core region, and this difference is mostly due to the effect of the scrape-off-layer but is also enhanced by grid resolution. Like the particle fluxes, the radial heat (conductive) fluxes are strongly suppressed towards the simulation boundaries in the annulus simulation. The level of the ion heat flux in the middle of the simulation domain is, however, equal in both simulations, whereas the electron heat flux at mid-radius is significantly lower in the annulus simulation due to the effect of the boundary conditions and grid resolution. Although the fluxes are larger in the full-torus simulation, the divergence of the electron and hydrogen radial energy flux is rather similar between the simulations, explaining the similar

evolution of the temperature profiles. The impact of grid resolution on turbulence dynamics is illustrated in Fig. 8, which shows the spatial scale of the electron density fluctuations to be much smaller in the higher-resolution full-torus simulation than in the lower-resolution annulus simulation. In conclusion, higher grid resolution and natural boundary conditions play an important role in improving plasma turbulence simulations. 6. Summary and conclusions The full f 5D gyrokinetic code ELMFIRE has been extended from a radial annulus to full radius version by including the magnetic axis and scrape-off-layer regions to the simulation domain. The scrape-off-layer plasma is bounded in the code by material surfaces along a radial wall and toroidal limiter, the latter of which is modeled as a toroidal chain of poloidal limiter plates for numerical reasons. In the implementation of the material boundary, the ion gyro-rings are let to penetrate into the radial wall, and the effect of the potential boundary condition in the gyrokinetic Poisson equation is enhanced by artificial manipulation of the ion polarization drift close to the wall and limiter plates. These two techniques, together with artificial suppression of the ⟨E⟩ × B drift next to the limiter plates, keep the simulations numerically stable with least effect on plasma edge structures. No other artificial tools, such as reflective boundary condition, decorrelation of particle motion and smoothening of the poloidal and toroidal plasma variations of the potential,

T. Korpilo et al. / Computer Physics Communications 203 (2016) 128–137

are needed in full-torus calculations. The applied manipulations, however, are not unambiguous at the wall boundary: the more aggressive the manipulation of the polarization drift, the lower the potential next to the wall surface, which improves the stability, but also increases the particle losses to the wall. The plasma center, in turn, is included in the simulation domain by using the magnetic coordinates which are singular at the magnetic axis. The effect of the coordinate singularity on the particle motion and the electrostatic field is avoided in ELMFIRE with the presented methods. From the numerical point of view, the extension of the simulation domain to the magnetic axis removes the inner radial boundary or one boundary prone to instabilities from the calculation. The predictive power of full-torus ELMFIRE simulations was examined by testing the code’s ability to reproduce experimental plasma equilibrium for the small FT-2 limiter tokamak. A FT-2 discharge was adopted as the basis of the simulation–experiment comparison because the configuration and parameters of this tokamak are well suited for ELMFIRE — only the experimental setup of two poloidal ring limiters was significantly different in the simulation that used the toroidal limiter configuration due to the numerical problems with the former. A direct comparison of the simulation results with the FT-2 data showed that the experimental density profiles were well maintained in the simulation and that the plasma current profile based on ASTRA predictions was correctly reproduced by the loop voltage close to the experimentally observed value. A poor power balance between ohmic heating, radiation cooling and radial energy transport losses was, however, observed in the simulation, which resulted in a steady increase of the plasma energy content or steady electron heating across the plasma. Simultaneously, the ion temperature profiles were modified substantially due to profile relaxation, further underlining the fact that the experimental equilibrium was not a numerical steady-state solution. A poor power balance and a strong relaxation of ion temperature profiles were also found in the simulation that adopted the initial parameters from another, lower-density FT-2 discharge. The result was verified by crosscomparison of the present code version and a previous one. These facts suggest that some plasma processes important for the energy losses in FT-2 are not included or modeled accurately enough in the present simulation model. This issue will be investigated further in the future. The cross-comparison of the code versions revealed also the impact of grid resolution and boundary conditions on turbulent transport. In other words, radial particle flux levels and radial electron thermal transport were significantly enhanced by highresolution grid and natural boundary conditions implemented in the full-torus simulation. Hence, the turbulence prediction capability of the present ELMFIRE version is clearly improved

137

over the previous versions, although it does not reproduce the experimental FT-2 profile equilibrium. The ELMFIRE results of scrape-off-layer turbulence simulations will be addressed in more detail in a future publication. Acknowledgments This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 under grant agreement number 633053 [and from Tekes — the Finnish Funding Agency for Innovation under the FinnFusion Consortium]. The views and opinions expressed herein do not necessarily reflect those of the European Commission. The work has been supported by the Academy of Finland grant 278487, RFBR grants 13-02-00614 and 15-02-03766, and High Level Support Team. CSC — IT Center for Science Ltd. and IFERC-CSC are acknowledged for allocation of computational resources for this work. Part of the results of this research has been achieved using the PRACE-3IP project (FP7 RI312763) resource Lindgren based in Sweden at PDC Center for High Performance Computing. References [1] S. Ku, C.S. Chang, P.H. Diamond, Nucl. Fusion 49 (2009) 115021. [2] P. Tamain, Ph. Ghendrih, E. Tsitrone, V. Grandgirard, X. Garbet, Y. Sarazin, E. Serre, G. Ciraolo, G. Chiavassa, J. Comput. Phys. 229 (2010) 361. [3] J.A. Heikkinen, S.J. Janhunen, T.P. Kiviniemi, F. Ogando, J. Comput. Phys. 227 (2008) 5582. [4] S. Leerink, V.V. Bulanin, A.D. Gurchenko, E.Z. Gusakov, J.A. Heikkinen, et al., Phys. Rev. Lett. 109 (2012) 165001. [5] J.A. Heikkinen, M. Nora, Phys. Plasmas 18 (2011) 022310. [6] T. Korpilo, J.A. Heikkinen, S.J. Janhunen, T.P. Kiviniemi, S. Leerink, F. Ogando, J. Comput. Phys. 239 (2013) 22. [7] J.A. Heikkinen, T. Korpilo, S.J. Janhunen, T.P. Kiviniemi, S. Leerink, F. Ogando, Comput. Phys. Comm. 183 (2012) 1719. [8] T. Takizuka, H. Abe, J. Comput. Phys. 25 (1977) 205. [9] S. Ma, R.D. Sydora, J.M. Dawson, Comput. Phys. Comm. 77 (1993) 190. [10] T.P. Kiviniemi, S. Leerink, P. Niskala, J.A. Heikkinen, T. Korpilo, S. Janhunen, Plasma Phys. Contol. Fusion 56 (2014) 075009. [11] D.V. Kouprienko, A.B. Altukhov, A.D. Gurchenko, E.Z. Gusakov, M.Yu. Kantor, S.I. Lashkul, L.A. Esipov, Plasma Phys. Rep. 36 (2010) 371. [12] G.V. Pereverzev, P.N. Yushmanov, IPP-Report, 2002, IPP 5/98. [13] M.Yu. Kantor, L.A. Esipov, D.V. Kouprienko, I.C. Nascimento, J.H.F. Severo, S.P. Yaroshevich, Proceedings of the 13th International Symposium on LASERAIDED PLASMA DIAGNOSTICS, Japan, 2007. [14] W.A. Houlberg, K.C. Shaing, S.P. Hirshman, M.C. Zarnstorff, Phys. Plasmas 4 (1997) 3230. [15] T.P. Kiviniemi, J.A. Heikkinen, A.G. Peeters, Contrib. Plasma Phys. 42 (2002) 236. [16] M. Kotschenreuther, G. Rewoldt, W.M. Tang, Comput. Phys. Comm. 88 (1995) 128. [17] N.T. Howard, C. Holland, A.E. White, M. Greenwald, J. Candy, Phys. Plasmas 21 (2014) 112510. [18] C. Holland, A.E. White, G.R. McKee, M.W. Shafer, J. Candy, R.E. Waltz, L. Schmitz, G.R. Tynan, Phys. Plasmas 16 (2009) 052301.