Developments of the culham 3D nonlinear MHD code

Computer Physics Communications 24 (1981) 437—439 North-Holland Publishing Company

437

DEVELOPMENTS OF THE CULHAM 3D NONLINEAR MHD CODE A. SYKES Cuiham Laboratory, Abingdon, Oxon, 0X14 3DB, UK (Euratom/UKAEA Fusion Association)

The structure of the 3D plasma simulation code NONLIN is outlined, together with a review of previous results. A new toroidal version, TORUS, is being used to study the behaviour of the internal sawtooth oscillation at high beta, and preliminary results are given. The conversion of NONLIN to run on the ICL Distributed Array Processor (DAP) is described, and it is shown that a very high performance is achieved.

1. Introduction The Cuiham 3D plasma simulation code NONLIN has so far been applied to straight systems in cases where the gross physics was believed to be insensitive to toroidal effects. Published work includes an explanation of the internal sawtooth in tokamaks [1] field reversal in pinches [2] and (in 2D helical form) studies of disruptions in tokamaks [3]. However, recent experimental work on ISX.B has produced high beta-poloidal plasmas with large toroidal shifts and unexpectedly good stability. Also, the JET apparatus has a very tight aspect ratio where significant toroidal effects can be expected. For these reasons the Cuiham code has been extended to toroidal geometry. Details of the code will be given, together with some results obtained from preliminary studies of the behaviour of the internal sawtooth oscillation for various values of beta-poloidal. NONLIN has also been programmed for a parallel computer of the array processor type (the ICL Dis~ tributed Array Processor, DAP). It is shown that this type of code is well suited to the DAP and a very impressive performance is obtained,

—

(p +~B2)b11} + V (pvpv1),

~-

—

~ip =

—

V• (pu) + V~(cVp), —

V (pu)

(2)

—

(‘y

l)p V

u

+~‘y—l~H,~r~j —f(xp

=

—

V XE,

E = —v X B

+

~

/= ~ X B

3. Applications at Culham

a

It is ambitious (and often unnecessary) to solve these full equations. The most common simplifica-

—pu1v~+ B1B1

0010.4655/8 l/0000—0000/$02.75 © 1981 North-Holland

(5

Our initial venture into plasma simulation was inspired by ORNL work [4], although we have always chosen to use the two-step Lax—Wendroff scheme. Since the first two equations are in conservative form, the 3 components of momentum and the density are conserved; however, the total energy is not, being determined by the balance of voltage input and applied losses such as radiation loss. For conducting shell cases, eq. (4) implies that the total B5 flux is conserved. The resistivity, viscosity and diffusion terms are solved in a single timestep dt, care being taken to maintain the appropriate conserved quantities, and to observe the restrictions on dr imposed by the various numerical stability conditions.

The MHD equations are solved in the following form:

-sa- —

3 (4)

2. Model

~—(pv1) =

(1)

438

A. Sykes / The Cuiham 3D nonlinear MHD code

tions are either to remove the inertia terms and/or to assume a given helicity, e.g. to just consider m = 2 modes. For each application we must decide which physical processes are important, treating these as accurately as possible; which are irrelevant, ensuring that computational time is not wasted; and if any approximations can be made to ease the computations. It is interesting to examine the various applications of the Culham NONLIN code in this light, a) Study of sawtooth phenomena in tokamaks (1976). The mechanism of the relaxation oscillation was then uncertain, and a full 3D treatment was necessary. However, neither toroidal effects or edge effects were considered important and so a straight geometry, square cross-section, wall-on-plasma model was used. Since only a basic understanding of an experimentally robust phenomenon was required, the ratio s of resistive times to Alfvén times was decreased considerably to s = 1000. This simulation, using a 14 X 14 X 10 mesh, was very successful. The value s = 1000 was large enough to distinguish inertial effects (time ~l), hybrid effects (time ~30), and resistive effects (time ‘~a1000). b) Study of field reversal in pinches (1977). This has very similar requirements to the sawtooth work, except that the HBTX pinch has low magnetic Rey. nolds number, s; in fact a complete experimental shot can be simulated. The use of a full 3D code was important for this problem; after the initial n = 1 helix which produces reversal, an n = 2 instability was seen to produce almost axisymmetric reversal, as is often observed experimentally, c) Study of disruptions in tokamaks (1980). Here a ‘wall-on-plasma’ model will not do: edge effects are known to be experimentally important. Also the onset of disruption is known to be sensitive to the size and behaviour of the m = 2 island: to represent this will require both a considerably finer mesh, and also a realistic thermal conductivity, with K11 ~ K1. A full 3D treatment to this accuracy was impractical on the available computers, and so it was decided to use the 2D form produced by assuming m = 2 helicity. The 2D equations were solved in a quadrant with circular section and a vacuum region and limiter. High parallel conductivity was achieved by integrating along the field line at each meshpoint, and setting the temperature at that meshpoint to a weighted average. In this way no detriment to K1 occurs.

This work gave good agreement with experiment: showing development of the m = 2 tearing mode until interaction with the cold limiter (or cold gas puffing region) leads to rapid growth of the island. At this stage other modes would have been excited; an important feature was that the central q was rapidly lowered. Representing the m = 1 sawtooth by an internal redistribution, the simulation showed that the m = 2 and m = I modes alone could lead to total disruption.

4. The effects of high-beta on the ‘sawtooth’ oscillation Recent experimental results from ISX-B [5] mdicate that for Ohmic discharges typical sawtooth behaviour occurs; but as neutral injection power is increased the period and amplitude of the X-ray signal first become enhanced, until at higher power the sawtooth can disappear, leaving prolonged precursor oscillations. Computational studies by Holmes et al. [6] solve a reduced set of equations and show that at high enough beta the shift of the current profile relative to the q = 1 surface stabilises the current driven m = 1, n = tearing mode island. However, their model does not include the pressure driven internal kink, which would presumably still be unstable for qaxis < 1. The Cuiham code TORUS (an adaptation of NONLIN to toroidal geometry) has been applied to this regime, using a 16 X 16 X 11 mesh on a 3 : aspect ratio torus of square cross-section with conducting walls. Runs have been made at values of beta

II 1

VE

I

He(~ty

(~3>S

1.6 ~I.

I I i

/

i

—

~

,,

— —

(p3)34’/.

—

I

100

200

300

400

Fig. 1. —ye helicity vs. time, for low and high values of beta.

A. Sykes / The Culham 3D nonlinear MHD code

corresponding to the two different experimental regimes, and preliminary results are presented here, The onset of the MHD instability is most clearly seen in the change in —ye helicity (V curl v), and fig. 1 shows the plots of —ye helicity, normalised to the kinetic energy at t = 300, against time for both a low and a high beta case. Both runs have the same total current, giving q(edge) = 2.7. The magnetic Reynolds number is typically S = 1000; Spitzer conductivity is used. The values of beta-f at t = 300 are 0.57 and 0.96, respectively. It is seen that the low-beta case performs regular sawteeth at a hybrid frequency whilst the higher beta case undergoes large, slow instabilities which evolve into a steady-state convective cell. The crudeness of the mesh, and the very low value of s chosen, prohibit any firm conclusions from being drawn; however, the suggested explanation is worthy of consideration namely that a) at low beta the m = 1, n = 1 instability quickly quenches itself by ensuring that q > 1 everywhere, and only reappears when the central q again falls below 1 due to plasma heating; b) at high beta the q < 1 region is not removed, and the instability remains in a mild nonlinear form, the associated convective cell removing any excess pressure. In this state the pressure is pushed into a horseshoe formation which, if rotated at the diamagnetic frequency, would give the appearance of ‘precursor’ oscillations.

439

written in simple structured form. As such it is ideally suited for conversion to DAP Fortran, which is a straightforward but powerful extension of ordinary Fortran. Although nearly all of the computer time used by NONLIN is spent in the utility routines AyE, Dlv, GRAD and CURL, it was necessary to convert the main timestep routine STEPON as well, since each entry to DAP involves conversion overheads. The gain in speed depends largely on how much of the DAP 64 X 64 processing array can be used. For example, a 16 X 16 X 16 3D problem fits exactly onto a single 64 X 64 DAP plane and so maximum efficiency is obtained: the speed gain is approximately 100 over the standard code running on the ICL 2976 with Optimising Compiler. Similar gains would be obtained for 21 X 21 X 18 or 32 X 32 X 16 meshes.

—

Acknowledgements The plasma simulation work has been done in collaboration with J.A. Wesson. The original conversion of NONLIN to DAP Fortran was done by Dan Oldfield of ICL.

References lii

5. Conversion to DAP

[2] [3]

The ICL Distributed Array Processor (DAP) is a 64 X 64 array of processing elements which perform computations simultaneously. The DAP itself replaces a storage unit on a 2900 series mainframe and is accessed by a program running on the host mainframe. NONLIN is an Eulerian, square mesh, explicit code

[4] [5] [6]

A. Sykes and iA. Wesson, Phys. Rev. Lett. 37 (1976) 140. A. Sykes and J.A. Wesson, Proc. VIII Eur. Conf., Prague (1977) p. 80. A. Sykes and J.A. Wesson, Phys. Rev. Lett. 44 (1980) 1215. J. Woolen, H.R. Hicks, G. Bateman and R.A. Dory, ORNL TM4784 (1974). V. Pare et al., paper 7P12, APS Conf., San Diego (1980). J.A. Holmes et a!., paper 5Q7, APS Conf., San Diego (1980).