P J Knight, C M Roach, A Thyagaraja, D J Applegate and N Joiner, UKAEA Fusion, Culham Science Centre

HPCx is the name of the UK's National High Performance Computing Service. It is a large IBM POWER5 cluster whose configuration is specifically designed for high-availability capability computing. The Engineering and Physical Sciences Research Council (EPSRC) is overseeing the project, on behalf of the UK scientific community. HPCx is a joint venture between EPCC (Edinburgh Parallel Computing Centre) at the University of Edinburgh and the Daresbury Laboratory of the Council for the Central Laboratory for the Research Councils (CCLRC). IBM (UK) Ltd has been chosen as the hardware supplier for the six-year duration of the project.

The late Nobel Laureate, Hans Bethe demonstrated that thermonuclear fusion powers the stars. In July of this year Cadarache in France was chosen to be the site for ITER (http://www.iter.org/), the world's largest machine dedicated to research into the production of energy from controlled thermonuclear fusion. ITER is designed to demonstrate that safe, clean electricity can be produced economically. When ITER begins to operate in about ten years it will aim to demonstrate that fusion can fulfil the energy requirements of the modern world without concomitant greenhouse gas production or need for disposal of long-lived radioactive wastes.

Achieving controlled thermonuclear fusion of light elements on Earth is not easy. In ITER and machines like it, a hot, ionised gas called a plasma must be confined for several seconds in a carefully designed 'magnetic bottle' well away from material walls. For fusion to take place in the deuterium-tritium fuel mixture (these are heavier isotopes of hydrogen, with one and two extra neutrons respectively), the temperature (>=100*10^6 degrees Celsius) must exceed by ten to twenty times that prevailing in the core of the Sun. The high temperature is needed to overcome the 'Coulomb barrier' to nuclear reactions between the deuterium and tritium nuclei. The magnetic fields (~5 Tesla) must be large enough to contain the plasma at a few atmospheres so that sufficiently many reactions can take place and the total power produced by them is significantly larger (>10 times) than that needed to heat the plasma. This plasma confinement approach leads to a doughnut-shaped device known by its Russian acronym, tokamak (toroidal, axisymmetric magnetic chamber).

The United Kingdom Atomic Energy Authority (UKAEA) carries out fusion research for the UK Government and the European Atomic Energy Community (EURATOM) at the Culham Science Centre in Abingdon, where research in experimental and theoretical physics and engineering tasks in support of ITER is undertaken. Our group's effort is devoted to the problem of understanding electromagnetic turbulence in plasmas, and its undesirable effects on enhancing energy losses (called 'anomalous transport') of the plasma far above those due to collisional processes.

Global tokamak turbulence calculations present truly 'grand challenges' to the most powerful computers in the world such as HPCx and the Earth Simulator in Japan. A whole range of linear and non-linear instabilities are involved. The problems are similar in complexity to those encountered in geophysical fluid dynamics and climatology, both in the enormous range of length and timescales involved and the number of dynamical degrees of freedom modelled. We are, in effect, attempting to 'arithmetize' plasma climatology in a tokamak with a supercomputer! We employ two complementary approaches to plasma turbulence modelling at Culham. The kinetic approach focuses on the microscopic spatial scale, exemplified by the Larmor gyration radius of charged particles in the magnetic field, whilst the 'continuum/fluid' method deals with more global scales right up to machine size and confinement time-scale, much longer than the turbulence timescale.

GS2 - a gyrokinetic code

Kinetic theories model the plasma as a collection of charged particles (~10^20) moving in response to self-consistently generated electromagnetic fields. Kinetic equations and Maxwell's equations are solved self-consistently to obtain the time evolution of the particle distribution functions in 6D (real space and velocity space) for each plasma species. Fortunately, in strongly magnetised tokamak plasmas it is possible to average the kinetic equation over the extremely rapid Larmor orbit motions to yield the 5D 'gyrokinetic equation'. This equation faithfully describes short perpendicular wavelength plasma turbulence, which is strongly suspected to underlie anomalous transport in magnetised plasmas.

GS2 [1] is a mature leading-edge plasma turbulence code, developed in the United States using F90 and MPI, to solve the non-linear gyrokinetic equation for each plasma species. GS2 is being used to study the microstability properties of spherical tokamak (ST) plasmas, including plasmas from the Culham MAST (Mega-Amp Spherical Tokamak) experiment. ST geometry presents a challenge for plasma theory, and has important influences on microinstabilities. GS2's domain is a 'flux-tube' sub-region of the tokamak plasma. Domain decomposition is performed in 5D, with care to minimise communications in the directions of the fastest processes along the magnetic field. GS2 has been widely exploited on supercomputers (including HPCx), and with judicious choice of the calculation grids scales efficiently with large numbers of processors.

Linear gyrokinetic calculations suggest the existence of key microinstabilities which may ultimately be responsible for anomalous transport of heat and particles. Magnetic field perturbations are found to be particularly important in the ST, and some instabilities are found to tear the equilibrium magnetic field [2]. Fully non-linear calculations are required to predict turbulence saturation levels, and to make contact with experimental observations.

Non-linear GS2 simulations on HPCx have calculated the saturated state of electron temperature gradient (ETG) driven drift wave turbulence in MAST plasmas. These calculations required 256 processors running for approximately 8 hours on HPCx, and predict a level of electron heat transport which is, remarkably, comparable to that measured. This result is of considerable interest as ETG turbulence has previously been dismissed as unimportant! Our recent HPCx non-linear ETG simulations for MAST [2], have revealed that radially extended 'streamer' structures in the electrostatic potential enhance electron heat transport. In the future it will be important to assess the robustness of these simulations to various approximations, and to find ways of 'tweaking' the plasma equilibrium to reduce the transport levels. We are also presently trying to understand the complicated physics mechanism that drives the instabilities that tear the equilibrium magnetic field lines. Even linear calculations of these modes are challenging computationally, and it will be important to make non-linear simulations to assess the level of transport arising from these modes.

CENTORI - a fluid code

Fluid models contain less detailed physics than gyrokinetics, although many observational phenomena on larger scales are within their ambit. These models represent a tokamak plasma using eight or nine 3D fields varying in space and time. Non-linear balance equations of particle number, momentum and energy for each species and a reduced set of Maxwell equations are derived by suitable averaging of the more exact kinetic description. CENTORI, the global fluid code being developed for use on HPCx, permits calculations of tokamak turbulence evolution and transport at reasonable resolutions (well beyond current experimental techniques but below those possible with gyrokinetics) for times approaching the typical confinement times.

CENTORI is designed to take full advantage of the parallel architecture of HPCx and Beowulf clusters. It is based on a highly successful serial Fortran 77 code [3] CUTIE (also developed at Culham by two of the authors), which described many experimental observations qualitatively with relatively low resolution calculations. CENTORI improves upon many of the approximations used in CUTIE, and employs a parallel implementation (MPI + FORTRAN 90) of the mixed semi-implicit pseudo-spectral/finite-difference method. Parallel supercomputing enables us to achieve much higher resolution than possible with CUTIE. CENTORI solves strongly coupled non-linear parabolic differential equations related to the well-known advection diffusion equation. A plasma-based coordinate system, based on the nested toroidal magnetic flux surfaces is used. The three spatial, plasma-based coordinates are: psi - labelling the nested flux surfaces, is a measure of the radial distance from the plasma centre to its edge; theta - the poloidal angle within a flux surface; and zeta - the toroidal angle along the torus. A predictor-corrector, semi-implicit finite difference scheme in psi is used with fast Fourier transforms in the two angular directions.

Early runs of CENTORI on HPCx have produced promising physical insights. Some typical outputs from the code, relating to UKAEA Fusion's own MAST experiment, already mentioned, demonstrate qualitatively the suppression of turbulence by the presence of a local bulk plasma flow in the toroidal direction. (The phenomenon is known as 'transport barrier formation' and is observed in many tokamaks.) Results were collected from another simulation, where CENTORI was used to model Alfvén wave propagation (so called after their Swedish discoverer Alfvén ) in the vicinity of saddle points ('X-point nulls') in the poloidal magnetic field. In this simulation (a two-dimensional, linear problem) we used 992 processors on HPCx with unprecedented resolution, and the results were in agreement with non-trivial exact solutions of the MHD wave equations [4]. Runs of this kind (only possible on HPCx!) have provided understanding of wave phenomena involved in particle acceleration and the possible loss mechanisms in X-point geometries.

The parallel solvers in CENTORI are being optimised with valuable code development support from Joachim Hein and Lorna Smith at EPCC as part of our EPSRC-HPCx contract. When fully developed and tested, CENTORI will be a powerful capability computing code. Global simulations with it on HPCx will provide new insights into both tokamak physics and plasma astrophysics.

Conclusions

A judicious combination of GS2 and CENTORI on HPCx and its successors will contribute to new understanding of existing tokamaks like JET (Joint European Torus), MAST and ultimately ITER. Theory and computation are expected to play an essential role in support of experiment and observation in finding a viable new solution to the world's energy problems and in understanding plasma phenomena in the cosmos.

Acknowledgements

This work was funded jointly by the United Kingdom Engineering and Physical Sciences Research Council and by EURATOM. The HPCx computer time was provided under EPSRC grant GR/ S43559/01. We would like to thank Bill Dorland for supplying GS2 and Joachim Hein and co-workers at EPCC for assisting in parallelization of CENTORI.

References

[1] http://gs2.sourceforge.net/
[2] C.M. Roach et al, to appear in Plasma Physics and Controlled Fusion (2005).
[3] M.R. de Baar et al. Phys. Rev. Letts, 94, 035002 (2005).
[4] K.G. McClements et al, 31st European Physical Society Conf. Proceedings, London, (2005).

This article originally appeared in the Autumn 2005 issue of Capability Computing, the newsletter for the HPCx community.  For more information about HPCx, visit http://www.hpcx.ac.uk/.