FUSION SIMULATION PROJECT (FSP)
Integrated Simulation & Optimization of Fusion Systems
Appendix. Overview of Fusion Science for the FSPI.

Introduction – Fusion Science Insights and Challenges
A. Insights from numerical simulation
B. Simulation challenges

II.

Fusion Science Basic Concepts
A. Basic equations
B. Reduced descriptions

III.

Microscopic Simulation
A. Introduction and summary
B. Neoclassical transport
C. Turbulence and anomalous transport
D. Non-local phenomena
E. Challenges for turbulence simulation

IV.

Macroscopic Simulation
A. Introduction and summary
B. Equations
C. Status
D. Challenges

V.

Plasma Edge Physics and Plasma-Wall Interaction
A. Introduction
B. Properties of the edge plasma
C. Computational challenges and codes for the edge region
D. Whole Device ModelingVI.

Equilibrium Evolution and Confinement
A. Introduction
B. Transport equations
C. Uses of transport codes
D. Challenges of transport modeling

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
1

VII.

External Sources
A. Introduction
B. Neutral beam injection
C. Edge particle sources
D. Radio-frequency wave heating and current drive
E. Summary: the role of source models in comprehensive integrated simulation

VIII.

Fusion Simulation Capabilities Status

IX.

Focused Integration Initiatives – Challenges and Opportunities
A. Plasma edge
B. Turbulence on transport timescales
C. Global stability
D. Whole device modeling

X.

Validation Requirements in Fusion Topical Areas
A. Transport
B. MHD
C. Radio-frequency heating and current drive
D. Edge/scrape-off-layer/divertor

XI.

Background – Orbits and Instability Mechanisms

XII.

Glossary

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP2

I. INTRODUCTION – FUSION SCIENCE INSIGHTS AND CHALLENGES
T his appendix summarizes the status of fusion science, emphasizing its theoretical andcomputational aspects. We wish to demonstrate first the scope and complexity of the physics andmathematics, second the progress that has been made up to now, and third the great challenges
and opportunities before us. We discuss the several areas of fusion physics and computation thatare referred to in the main report, beginning with a brief introduction to the basic equations andconcepts in Sec. II. We follow this with more detailed surveys of the main areas of concern in
Secs. III-VII. In Sec. VIII we summarize the 50 or so principal codes in use today. In Secs. IX
and X we explore the focused integration initiative examples and validation requirements that are
introduced in the main body of the report. In Sec. X we give a very brief review of someelementary plasma physics concepts. Finally, in Sec. XI, there is a glossary. Words introduced initalics are defined there.
Before going into the various topics, let us begin by examining our ultimate goal from the pointof view of past accomplishments in computational plasma physics and some challenges that we
expect will lead to future insights.
Numerical simulation and modeling, at varying stages of integration, have contributed greatly toinsights regarding the behavior of magnetically confined plasmas. This has resulted from thecontinual interaction among computation, theory, and experiment. As the Fusion Simulation
Project (FSP) develops, and wider ranges of physical processes are integrated into thesimulations, there will be greater and greater reliance on massive computations. At the sametime, as pointed out in the main text, the theoretical interpretation and experimental validation
efforts will continue to be crucial to the overall success of the project.
We present a few examples to illustrate how this process works and how computation has
contributed and will contribute to our insight and understanding.
A. INSIGHTS FROM NUMERICAL SIMULATION
1. Global stability and major disruptions in tokamaks
A major disruption is the name given to an abrupt abnormal termination of a tokamak discharge
with total loss of plasma current and confinement. This is thought by many to be the major
obstacle that could potentially upset the development path of the tokamak as a practical energy
source. Besides terminating the discharge, the disruption has the adverse effects of producing
large thermal loads on the divertor and large electromagnetic forces on the vessel wall. (Seecutaway view in Fig. III.1 of the main report.) Disruptions also have the potential of accelerating
electrons to multi-MeV energies, and these can cause additional severe damage to the first wall.
Computer simulations have been crucial in identifying specific sequences of events that can cause
disruptions. An excellent example of this are the simulations that explained the major disruption
in the record-setting TFTR discharge that produced 10 MW of fusion power. [E. Fredrickson,
W. Park, et al, “High-beta disruption in tokamaks”, Phys. Rev. Lett, 75 1763 (1995)]. It was
shown that under the conditions of that discharge, it was energetically favorable for the plasma

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
3interior to deform into a long-wavelength helical structure. In the local regions of the torus
where the unfavorable curvature of the helical structure aligned with the curvature of the torus
itself, the plasma pressure could very easily cause instability. It was the interaction of thislocalized ballooning instability with the interior long-wavelength helical structure that destroyed
the nested magnetic surfaces and precipitated the disruption. Figure A1 depicts a simulation ofthis process, with perturbed field lines and poloidal cross sections highlighted in colors
representing mode amplitude.

Fig. A1. An internal helical structure interacting with the toroidal geometry causes localized instability that can lead to a major
disruption. [Courtesy Wonchull Park, PPPL]

There are other causes of the major disruption, which we are just beginning to understand andare finding ways to control. Integrated simulations of the kind planned in this project will be
instrumental in identifying and understanding these mechanisms. (See Sec. IV.)
2. Design of Three Dimensional Toroidal Configurations
Another example with already some degree of integration is the design of three-dimensional (3D) stellarator configurations. Stellarators are advanced toroidal magnetic confinement devices
that utilize 3-D plasma shaping to achieve steady-state operation and to optimize plasma
confinement and stability. Among these, the design of compact (or low aspect-ratio) stellarators
is one of the new, challenging research areas of interest in magnetic confinement. The challenge
is that maintaining good transport and stability properties at low aspect ratio is generally moredifficult than at high aspect ratio.
However, careful control of the 3-D shape appears to offer the flexibility to achieve goodtransport and stability characteristics. This insight has been realized through the development ofsophisticated optimization codes that numerically explore 30-40 dimensional parameter spaces(the stellarator’s outer flux surface shape is described in terms of 30-40 Fourier modes) and arrive
at shapes that minimize a range of desired equilibrium, transport and stability physics targetfunctions [A. S. Ware et al., “High-beta equilibria of drift-optimized compact stellarators”, Phys.

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
4

Rev. Lett. 89, 125503 (2002)]. Typical stability targets are the Mercier criterion, resistive
interchange, high-n ballooning, and kink mode stability.

Fig. A2. Quasi-axisymmetric National Compact Stellarator Experiment (NCSX). [Courtesy G. H. Neilson, PPPL]

These optimization procedures have led to unexpectedly robust configurations such as the quasiomnigeneous or quasi-poloidal compact stellarator (QPS) being designed at Oak Ridge National
Laboratory and the quasi-axisymmetric National Compact Stellarator eXperiment (NCSX)
device [http://www.pppl.gov/ncsx/] currently in the design phase at Princeton Plasma Physics
Laboratory. Figure A2 is a plan view of NCSX, showing the outer closed flux surface and themain magnetic coils. The integration initiative could lead to the inclusion of more sophisticatedand more compute-intensive particle-transport, 3-D stability, and coil-reconstruction modules
directly in the optimization. This could yield yet more surprising and robust three-dimensional
configurations suitable for reactor-grade plasma confinement.
3. Turbulence Regulation by Self-Generated FlowsA major step in tokamak research was the experimental realization that the High Mode (HMode) of confinement in tokamaks was accompanied by a sudden increase in the primarily
poloidal rotation of the plasma in the region of steep pressure gradient at the plasma’s edge.
Almost contemporaneously, early computer calculations were showing that turbulence could
spontaneously generate large radial scale flows, which formed a transport barrier [A. Hasegawa
and M. Wakatani, “Self-organization of electrostatic turbulence in a cylindrical plasma”, Phys.
Rev. Lett. 59, 1581 (1987)]. This insight was later confirmed experimentally with the discovery
of internal barriers, confirmed computationally with the almost universal observation of
Reynolds-stress-generated flows in fluid, gyrofluid and gyrokinetic simulations, and put on soundtheoretical foundations by several thorough analytical studies.
FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
5

Fig. A3. Computational evidence of flow regulation of turbulence [Courtesy Z. Lin, PPPL/UC-Irvine]

In particular, the computations showed the effectiveness of these flows to regulate the turbulencein terms of scale size, fluctuation levels, diffusivities and therefore the transport. Recent
theoretical effort motivated by differences in the characteristics of the turbulence obtained fromfluid and kinetic models led to the identification of a residual flow damping mechanism, which
in turn gave rise to an improvement of the computational methods. Further theoreticaldevelopment has also led to the discovery of long-lived, fine radial scale flows nonlinearly
generated by the turbulence or zonal flows whose properties and consequences were confirmed in
the details by the computational models. Figure A3 represents a microturbulence simulation.
Poloidal cross sections at bottom left and right show flow patterns without and with zonal flows[Z. Lin et al., “Turbulent Transport Reduction by Zonal Flows: Massively Parallel Simulations,”
Science 281, 1835 (1998)].
We anticipate that these kinds of insights coming from the synergy among experiment, theory,and computation, exemplified by the role of flow in regulating the turbulence, will be one of the
distinguishing characteristics of the integration initiative. (See Secs. IV-VI.)
4. Full-wave modeling of radio frequency heated, multidimensional plasmas
A major challenge in understanding waves for fusion is that they can dramatically change
character after launch, i.e. undergo mode conversion. At a given frequency a uniform magnetizedplasma supports several different kinds of waves, with very different characters with respect topolarization, speed and wavelength, such as fast magnetosonic waves, slow ion-cyclotron waves
and ion Bernstein waves (IBW). (See Sec.VII.) In a non-uniform plasma these different waves can
FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
6

convert at isolated spatial locations, the classic example being conversion of the fastmagnetosonic wave to an IBW at the 2-ion hybrid resonance.
Understanding this process in the realistic 2- and 3-dimensional geometries of toroidal magnetic
confinement devices has required computational modeling and only recently has it become
possible to fully resolve the mode conversion layer, for example with the All Order Spectrum
Algorithm or AORSA code [E. F. Jaeger et al., “Advances in full-wave modeling of radio
frequency heated, multidimensional plasmas”, Phys. Plasmas 9, 1873 (2002)].

Fig. A4 Mode conversion to fast ion cyclotron wave revealed by AORSA computations. [Courtesy D. Batchelor, ORNL]

The surprise is that these realistic-geometry and well resolved calculations show that with full
poloidal field the fast magnetosonic wave converts dominantly to a slow ion-cyclotron wave, notan IBW. This had been hinted at by F. Perkins who did an approximate 1-D calculation in 1977
showing that both types of conversion could, in principle, occur [F. W. Perkins, “Heating
Tokamaks via the ion-cyclotron and ion-ion hybrid resonances”, Nuclear Fusion 17, 1197
(1977)]. But the 1-D model did not show that one could actually get a wave in a real 2-D
plasma, launched from a real antenna, to the conversion layer with the right conditions forconversion. Nor did it tell how important each process would be. The 2-D calculation gives thecomplete, quantitative picture, but the 1-D calculation gives the qualitative paradigm by which
one understands the complete results. Figure A4 shows wave amplitudes in a cross section of the
C-Mod tokamak, with details in the conversion layer expanded in the right-hand diagram.
As with the other examples cited, RF calculations with ever improving resolution and ever tighter
coupling to more comprehensive plasma models such as those contemplated as part of thisintegration initiative are expected to lead to even more surprises.
FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
7

B. SIMULATION CHALLENGES

1. Modeling edge physics
Overall transport and confinement are determined to a great extent by the height of thetemperature pedestal at the plasma edge. This temperature is itself governed by phenomena thatstart at the wall, and propagate through the separatrix region and into the edge region proper. It
is expected that coupled and complex models of particle and heat transport, neutral and impurity
fluxes, edge-gradient induced magnetohydrodynamic (MHD) instabilities, and turbulence in asingle computational edge framework will be required to pin down which of these mechanisms,
either by itself or in combination with another, regulates the pedestal height. (See Sec. IXA.)
Coupling to core plasma models will probably start with the integrated edge models providing
the pedestal temperature as a boundary condition for the core before more complete integration
becomes feasible through a combination of physics, algorithmic, and hardware advances.
2. Modeling turbulence on transport time scales
An integrated model of turbulence on transport time scales is a daunting physics and computationaltask. It is nevertheless deemed feasible at several levels, each exploiting separation of space and timescales appropriately. (See Sec. IXB.) One approach consists of iterative solutions of the macroscopictransport equations with occasional updates of the profiles and diffusivities transmitted to and fromthe microscopic turbulence calculations. Another entails direct coupling of transport equations togyrofluid-type microscopic models using a physics-based, yet reduced, description of the turbulence.
The result of high confidence integrated modeling of turbulence and transport might be thediscovery, through computations, of new favorable operating modes, such as a new H mode, withthe ultimate outcome being the determination of transport from first principles.
3. Coupling of electrons and ions in core turbulence
Integration will enable examination of the combined role of turbulence on space and time scales
associated with the electrons and ions. (See Sec. IIA.). To date realistic and tractable numericalsimulations of each exist separately. The expected advances in integration algorithms to bridge
disparate time and space scales will enable incorporation of the dynamics of each species withrealistic mass ratios, and for realistic scale sizes with respect to experiments, into one
computational model that covers the full range of relevant turbulence. The benefits of thesecomputations will be to determine the extent to which electron and ion turbulence overlap interms of spectral range and to what extent one range of scales influences and/or regulates thetransport associated with the other.
4. Extended MHD
Integration will facilitate extensions to MHD computations beyond the conventional ideal andresistive models. (See Sec. IV.) Routine incorporation of two-fluid effects will determine for
instance whether whistler wave physics included through the Hall term will play as important arole in phenomena leading to plasma disruptions as it does for magnetic reconnection in space
and astrophysical plasmas. Similarly, the routine inclusion of diamagnetic effects associated withFESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
8the non-scalar nature of the ion pressure tensor will reveal whether the pressure-induced plasma
rotation provides a way to control MHD activity which is as effective nonlinearly as it is linearlyfor realistic toroidal plasmas. Moreover, the inclusion of minority ion species with nonMaxwellian populations either as gyrofluids or full-fledged particle species will enable extended
MHD models to describe the essential macroscopic physics of burning plasmas. Such features
will provide new intuition regarding the effect of energetic particles on large-scale MHD
instabilities and vice-versa.
5. The burning plasma experiment
Perhaps the greatest advances afforded by integrated modeling will be realized for burning
plasmas. The grand challenge in the world fusion program is the burning plasma experiment, anecessary predecessor to a practical power reactor because, by its very nature, it presents a new
category of technical issues. All other fusion experiments obtain the required plasma
temperatures by applying energy sources from outside the plasma. A burning plasma – a selfsustaining energy source – maintains them through the action of the ongoing fusion reactions,which produce energetic alpha particles and high-energy neutrons. The charged alphas aretrapped by the magnetic field and transfer their kinetic energy to the plasma, thus maintaining
the temperatures necessary for fusion to continue.
With self-heating as the dominant plasma heating mechanism, new plasma processes and effects
will arise. The high flux of energetic particles will impact the plasma wave structure, alter theplasma pressure and current profiles, and produce, through alpha-particle and neutron collisions,a rich source of wall interaction phenomena. Most seriously, all these effects will be strongly
coupled and must be understood and managed in an integrated fashion to insure the stability
and success of the experiment. The fully integrated fusion simulator will be targeted to model
these processes and their consequences, thereby providing the essential insights to guideexperimental programs, optimize machine design, and deepen our understanding of thefundamental science.

II. FUSION SCIENCE BASIC CONCEPTS
The fundamental problem in fusion science is to understand and control the physical processes
that determine the balance between heating and fueling, on the one hand, and loss ofconfinement, on the other. The losses can come from macroscopic processes leading to breakupof the nested magnetic flux surfaces and disruptions or from microscopic processes, leading totransport across the flux surfaces. Figure 5A shows a tokamak cross-section, schematicallyillustrating some of the processes that will be discussed in the following sections. A cutaway viewof this device is shown in Fig. III.1 of the main report.

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
9

ECRF
launcher

ICRF or Lower
hybrid launcher

Resonance layer ω = ΩFig. A5. Tokamak cross-section with launchers, a schematic wave and resonance layer. Here Ω = ω c . Also shown is a chain ofmagnetic islands.

We begin in Part A by introducing the basic equations and concepts of plasma physics, which are
the foundation of the subsequent sections, and which give rise to the time-scale chart appearing
as Fig. III.2 of the main report. Reproducing this chart in Part B, we briefly outline these timescales.
A. BASIC EQUATIONS
The time-dependent processes governing the stability and confinement of the plasma aregoverned by the Maxwell-Boltzmann system of equations, which we write in MKS units. We
begin with the distribution function f a ( x,v,t ) , which is the density of charged particles of
species a in a six-dimensional phase space. This function evolves in time, t, according to theBoltzmann equation,

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
10

∂f a
q
+ v ⋅ ∇f a + a [E + v × B] ⋅ ∇ v f a = C ( f a ) .
∂t
ma

(A-1)

The characteristics of the left-hand side are the particle orbits, the second term, v ⋅ ∇f a , being
the convection in space and the third term, (qa ma )(E + v × B) ⋅ ∇ v f a , being the acceleration invelocity space arising from the Lorentz force, where E and B are the electric and magnetic fields,
and qa and ma are the charge and mass of species a. Solution of the equation can begin withintegration along the characteristics. (A brief Introduction to these is given in Sec. XI.) This iscomplicated, however, both by the right-hand side of the equation, which represents all sources,
sinks, and collisional dissipation, and by the nonlinearity of the whole system. That is, the fields
themselves include not only the externally generated B, but also spontaneous perturbations andinjected waves that depend on the particle distributions through the Maxwell equations,∇ × B = µ0 J + µ0ε0
∇×E =−

∂E
,
∂t

∂B
,
∂t

∇ ⋅ B = 0,

(A-2)

∇ ⋅ E = ρ ε0 ,(A-3)

where J and ρ are the sums of individual species current and charge densities,

J( x,t ) = ∑ J a ( x,t ) = ∑ qa ∫ d 3 vvf a ( x,v,t ) ,(A-4)

ρ( x,t ) = ∑ ρ a ( x,t ) = ∑ qa ∫ d 3 vf a ( x,v,t ) .

(A-5)

a

a

a

a

Theoretical and computational challenges. Equations (A-1)-(A-3) generate an enormous array ofphenomena. We shall touch upon the more important of these in the following sections. The
sources of this complexity are:
1. The complex structure of the confining magnetic fields, which give rise to multiple
classes of orbits even in the absence of perturbations,2. The long-range nature of the electromagnetic forces, which give rise to collective
phenomena, both externally driven and spontaneously occurring waves and instabilities.
3. The great range of space and time scales and extreme anisotropy of both individual
particle motions and plasma waves, extending from the size of the device (a few meters)
and the duration of the discharge (exceeding several minutes in a recent record-setting
experiment) all the way down to the electron gyroradius and gyrofrequency, and
4. Dissipation, i.e., collisions among ion species and electrons, as well as those with neutrals(sources and sinks), all buried in the term C ( f a ) on the right-hand side of Eq. (A-5).
Even when the collision term itself is negligibly small, dissipation occurs and plays amajor role, owing to wave-particle resonances, and nonlinear decorrelation of waves.

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
11

An electron in a uniform magnetic field B orbits the field line in a perpendicular plane (Part B,
below) with gyrofrequency ω ce = eB me ≈ 5 ×1011 s-1 a n d
withgyroradius−5ρ ce = v ω ce ≈ 6 ×10 m . (These numerical values are for a typical plasma temperature of 3 keV
and a typical magnetic field strength of B = 3 T.) Frequencies and wavelengths for the highest
range of radio-frequency wave heating are of the same order. (Ions gyrate with frequency lower
and radius larger by me mi and mi me , respectively.) The same electron moves freely parallelto the field, but in the curved and spatially varying fields of a torus, it is subject to additional
slower drifts and accelerations on the whole device scale, of order 2πR ≈ 10 m in present-day
machines, where R is the major radius.
Thus, our phenomena extend over a range of 9-10 orders of magnitude in time and,perpendicular to the magnetic field, 4-5 orders in space. For much of our work, we can also
average over the gyro-motion, reducing our equation to five phase-space dimensions (gyrokinetic
theory), and increasing the time step by another order of magnitude. Even so, one estimate of thecoverage needed for a full-discharge simulation would require 1011 phase-space points at 10 8time steps. Not only does this far exceed our present capability, but also we would not be able tostore or analyze the output from such a calculation.
Actual plasma computations, therefore, are carried out by means of ordering schemes leading toapproximations that focus on particular scales and are valid for particular phenomena. Beginningwith Eqs. (A-1)-(A-3), the equations to be solved are obtained by averaging over the finer space
and time scales and applying some prescription for closure. For example, in a toroidal device,
particles drift away from magnetic surfaces at a rate slower than the gyration by factors of orderρ* = ρ ci a , where a is the minor radius. Relatively infrequent collisions or decorrelations then2

allow transport losses described by a random walk with diffusivity of the form D ~ (δx ) δt
where δx and δt, the characteristic step size and collision time, exceed the shortest scales given
above. If D and related quantities can be simply characterized, then the resulting transport
equations are solved on still longer time and space scales. Whole-device simulation currently
proceeds in this manner, but the “if” in the preceding sentence is a big one.
Much progress has been made in the major subfields of plasma theory and computation, which
we describe in the following sections. Figure III.2 in the main report, repeated as Fig. A6 in thenext section, shows a chart of the principal frequency ranges and the phenomena associated withthem.

The great challenge and opportunity before us is to further develop and combine these in a tractableway, leading to whole-device simulation that contains complete and reliable models of all the relevant
physics.
B. REDUCED DESCRIPTIONS
The codes currently in use are each based on different approximations to the full set of plasma
equations (A-1) to (A-3) that isolate phenomena in a more restricted range of frequencies.
Together, they span the ranges of frequency described above. A summary of the four major codegroups presently in use and the timescales being addressed by them is given in Fig. A6.

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
12

Typical Time Scales in a next step experimentwith B = 10 T, R = 2 m, ne = 1014 cm-3, T = 10 keV
SAWTOOTH CRASHELECTRON TRANSITΩce-1 ωLH
Ωci-1
-10
-8
10
10
-1

TURBULENCE

τA

10-6

ISLAND GROWTH

10-4

Neglect displacement
Single frequency
current, average overand prescribed
plasma background gyroangle, (some)with electrons
RF Codeswave-heating andcurrent-drive

ENERGY CONFINEMENT

10-2

CURRENT DIFFUSION

100

102

104 SEC.

Gyrokinetics Codes

Neglect displacement
current, integrate overvelocity space, averageover surfaces, neglect
ion & electron inertia
Extended MHD Codes Transport Codes

turbulent transport

device scale stability

Neglect displacement
current, integrate overvelocity space, neglect
electron inertia

discharge time-scale

Fig A6. Summary of the four major code groups and the timescales being addressed.

The RF codes (external sources) address frequencies of order the ion cyclotron frequency, Ω ci,
and above, up to the electron cyclotron frequency Ωce. (This notation is used interchangeably
with ωci and ωce.) They assume a single, fixed frequency set by the oscillator, and solve for the
linear and quasi-linear response of a given background plasma to the imposed electromagnetic
fields at this frequency. They are used to calculate wave-heating and current-drive by radiofrequency (RF) sources.
The gyrokinetics codes (microscopic modeling) are based on an analytic averaging of the fullkinetic equation to eliminate the fast gyro-orbit frequencies from the equations. They and theother codes described below also neglect the displacement current term in Maxwell equation (A2) to remove light waves from the system. These reductions allow them typically to take time
steps about 10 times longer than the ion cyclotron frequency (whose motion is analyticallyaveraged over) and are normally run for 103 to 104 time steps to calculate stationary turbulent
fluctuation levels. Recent additions to these codes to include some electron timescale phenomenabring in the electron transit time, which lowers the maximum time step.
The Extended MHD codes (macroscopic modeling) are based on taking velocity moments (orweighted integrals over velocity space) of Eq. (A-1) to obtain fluid-like equations that describe
large-scale phenomena more efficiently. They aim to resolve phenomena occurring on the Alfvén

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
13

transit time, τA, although most codes are at least partially implicit to avoid a strict restriction on
the time step based on this. These codes can normally run 104 to 105 time steps to address MHD
phenomena such as sawteeth and island growth.
The 2-D transport codes (equilibrium evolution and confinement) further simplify the equations
by eliminating the Alfvén waves and by averaging plasma properties over 2-D magnetic surfaces.
These codes are very efficient and can take many long time steps to model the entire discharge.
However, the 2-D edge transport codes need to resolve the parallel dynamics and use a time-stepbased on the parallel sound-wave propagation near the edge region.
The following sections describe each of these code groups in more detail.

III. MICROSCOPIC SIMULATION
A. INTRODUCTION AND SUMMARY

This section deals with the direct calculation of the microscopic processes that affect the qualityof a magnetic field configuration as a thermal insulator. Loss of plasma particles and energy
across field lines results from at least three categories of microscopic phenomena: diffusion andconvection based on individual orbits and collisions (classical, neoclassical), anomalous diffusionand convection from turbulent microinstabilities (usually thought to be dominant), andphenomena that are instantaneous and non-local.
• Neoclassical codes for both 2- and 3-D configurations are available subject to certain
approximations. Additional work will be required to implement fully the existing
theories and to account for additional effects in the inter-operative environment.
• Micro-turbulence-driven anomalous transport will be dealt with on three levels: a) finescale stand-alone gyrokinetic simulations of core plasma turbulence will continue to be
developed, interpreted by theory, and compared with experiment to firmly validate thefundamental theory and benchmark the other descriptions; b) reduced simulations will
be coupled directly to the transport equation solvers. In some cases, it may be possible tocouple directly the full turbulence simulations with the transport equation solvers (atopic of one of our focused initiatives); and c) algebraic models of transport will continue
to be developed, again informed by theory, experiment, and the just-describedsimulations. When supported by turbulence simulations for selected cases, they will
provide the most rapid parameter scans• Work will continue on development of models of observed rapid and non-local
phenomena (e.g. avalanches, radiative transport, or global magnetic interactions) that do
not fit into the diffusive-convective approach. The architecture must include provision
for these from the outset, so they can be incorporated when available.
A measure of success for these efforts and their validation will be the ability to predict and modela transport barrier, a region of steep gradients where turbulence is suppressed. This will involve
integrated treatment of the radial electric field as also mentioned in Secs. V and VI.

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
14

B. NEOCLASSICAL TRANSPORTAt the most basic level, transport is a random-walk process, in which a particle drifts slowly awayfrom a magnetic flux surface and from time to time suffers a collision that transfers it to a new
2
orbit. The resulting diffusivities have the form D, χ ~ (δx ) δt where δ x and δ t are thecharacteristic step size and collision time. The detailed task of neo-classical theory is to computethe flux-surface averaged fluxes (see Sec. VI, Eqs. (A-24)-(A-25)) by solving a gyro-averaged formof the Boltzmann equation (A-1), called the drift-kinetic equation. Closure is provided in anordering in which the collision operator is dominant, and the step size (drift orbit) is small
compared to the plasma gradient scale length. (Classical transport, in which the step size is of theorder of the gyroradius ( δx ~ ρ a ), is small and most often neglected in fusion calculations.)
Neoclassical theory is highly developed in the limit that finite orbit-width effects (violation of theabove ordering) and energy scattering are unimportant. Codes for both 2- and 3-D
configurations are available subject to these approximations. Additional theoretical work andextensive code development will be necessary in order to account for a number of importantphenomena associated with neoclassical theory such as: neoclassical transport when barriers turnoff the anomalous transport, bootstrap effects (neoclassical currents driven by gradients), nonlocal orbit effects such as potato orbits (orbits near the plasma center), self-consistent electric field
effects in non-axisymmetric configurations, impurity transport etc. Additional work will be
required to fully implement the existing theories in the “interoperative” environment,
particularly in 3-D configurations (stellarators) and in tokamaks with transport barriers. To treat
the mutual interaction of islands or stochastic magnetic fields with neoclassical driven transportand currents will require extensive further development.
C. TURBULENCE AND ANOMALOUS TRANSPORT
Gradients of plasma density and temperature provide the free energy to drive micro-instabilities
that lead to anomalous transport. These waves typically propagate with phase velocities*
perpendicular to the field of order v *a = ρ a v a d ln n e dr (see Sec. XI, Eq. (A-30)), or vTa<< v awhere the temperature gradient replaces the density gradient. This gives rise to the characteristic**
diamagnetic frequencies ω *a = k⊥ v *a and ωTain a wave φ˜ ~ exp(ik ⊥ ⋅ x ) , where k⊥ ~ ρ −1
= k⊥ vTaa ,for a = i (ion waves) or a = e (electron waves). The waves are known generically as drift waves,particular examples of which are the ion temperature-gradient (ITG) modes and electron
temperature-gradient (ETG) modes.
To evaluate these we begin with the Maxwell-Boltzmann system, Eqs. (A-1)-(A-5), neglecting
displacement current and charge density (quasi-neutrality). Averaging over the gyro-angle, we
obtain the guiding center distribution Fa ( x,v,t ) = F0a ( ρ,W )(1+ qa φ˜ Ta ) + h˜ a ( x,W , µ,t ) in a 5-

dimensional phase space, where the velocity variables are the magnetic moment µ = mv ⊥2 2B
2and kinetic energy W = mv 2 2 = µB + mv|| 2 . It is convenient to separate out the background
distribution F0a ( ρ,W ) and the so-called adiabatic part F0a ( ρ,W ) qa φ˜ Ta . (The latter is actually
the lowest order solution in one important limit.) The remaining part of the distribution, h˜ a

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
15

obeys the gyrokinetic equation, in which the convective term in (A-1) is replaced by a similar termwith gyro-averaged particle drifts replacing the velocity v:∂h˜ a∂f ∂χ˜+ v χa + v da + v||bˆ ⋅ ∇h˜ a = −v χa ⋅ ∇f 0a − qa 0a∂t∂W ∂t+ collisions + sources/sinks,

(

)

(A-6)

where bˆ points in the direction of the equilibrium magnetic field, v da is the curvature- and gradB drift, and the E×B drift is combined with transport along perturbed magnetic fields lines andthe perturbed grad-B drift asv χa = bˆ × ∇χ a B ,

where


χ a =  φ˜ −v ||A˜||g+ µB˜|| / qa

(A-7)
 B.

g

(A-8)

Here, we represent the fluctuating fields in terms of the potentials (φ˜, A˜ || ) plus B˜ || and denote thegyro-average by the angle brackets … g . The principal nonlinearity of Eq. (A-6) is theperturbed convective term v ⋅ ∇h˜ .
χaa

To illustrate the method and make contact with theory, we note that the gyrokinetic ordering
requires a separation of scales that lends itself to an eikonal representation, i.e., a perpendicularexpansion of fluctuations as exp(iS ) where ∇S = k ⊥ may be interpreted as the perpendicularwave vector. In this case, Eq. (A-8) may be evaluated asχ˜ a = J 0 (ba )(φ˜ − v|| A˜ || ) + J1 (ba ) ba (µB˜ || /qa ) .

(A-9)

The Bessel functions J0 and J 1 have argument ba = k ⊥ v ⊥ ω ca . Of course, in a non-uniform
medium the perpendicular wave number is really a differential operator. (In order to evaluate
these functions efficiently, the codes we will describe make extensive use of numerical Fouriertransformation and its inverse.) In particular, the gyro-averaged wave potential is given by1
φ˜ a ( x ) =∫ dζ φ˜ (x + ρa ) = J0 (k⊥ ρa )φ˜ (x) ,2π(A-10)

where ζ is the gyro-angle, as illustrated in Fig. A7.
Finally, the Maxwell equations are recast in potential form, e.g.,∇ 2 A|| = −µ0 ∑ qa ∫ d 3vv|| J 0 (ba ) h˜ ( a) ,awhere the right-hand side is the parallel current.

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
16(A-11)

We draw attention to the fact that solving the continuum Eq. (A-6) is equivalent to following
gyro-averaged orbits of an ensemble of discrete particles,

Fa = ∑ w i ( t )δ ( x − x i ( t ))δ ( v − v i ( t )) ,(A-12)

i=1,N

with suitable particle smoothing, Monte Carlo treatment of collisions, sources and sinks.

∼
φ

ρa
ζ
x

Fig. A7 Gyro-orbit average of perturbed potential.

The system consisting of Eq. (A-6) and the Maxwell equations, Eqs. (A-2)-(A-3) (neglecting thedisplacement current,) drives turbulence on scales ranging from the electron to the ion
gyroradius, as noted at the beginning of this section. The linear and nonlinear growth anddamping mechanisms have been well delineated theoretically, and a number of codes are able tocompute the evolution and saturation of the waves and the resulting flux-surface-averaged transport
fluxes (see Sec. VI). The codes are of four types. Each may be either flux-tube (localized
perpendicular to the field) or global (covering a major fraction of the plasma radius). Also, eachmay be either continuum (solving Eq. (A-6) on a fixed Eulerian grid), or particle-in-cell (or PIC,
following an ensemble of gyro-averaged particles in a Lagrangian formulation) as in Eq. (A-12).
All four types of code have advantages and disadvantages.
The SciDAC numerical turbulence project supports a two-by-two matrix of codes:
Flux tube
GlobalContinuum
GS2
GYRO

PIC
SUMMIT
GTC

There are both continuum and PIC because continuum is currently the most developed (kinetic
electrons, important even for waves on the ion gyro-scale, and perturbed magnetic fields) while
PIC may be ultimately more efficient. A flux tube is described in coordinates that are extended
along field lines but are localized in perpendicular directions. There are high-resolution flux-tube
codes to support multiple space and time scales (electron and ion gyro-physics) and validation of
FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
17

local turbulence physics (rapid parameter scans). Global codes, which provide calculations over a
major fraction of the plasma radius, are needed to account for extended profile effects and aremost likely the best bet for coupling to the integrated simulation.
The plasma turbulence is quasi-2-dimensional because of the plasma anisotropy. The current
state-of-the art in global PIC simulations with GTC, without including full electron dynamics ormagnetic perturbations, is the following. The resolution requirement along B is determined by
the equilibrium structure. The structure across the field is determined by the microstructure of
the turbulence, ~ ρ ci for ion temperature-gradient (ITG) m o d e s . This requires2~ 64 × ( a ρ ci ) ~ 2 ×10 8 grid points and 8 particles per spatial grid point, where a is the minor
radius of the machine. This leads to ~1.6×109 particles and 1 terabyte of RAM for
600 bytes/particle. This resolution is currently achievable. The time scale is on the order of
a v i ~ 1 µsec for 10 time steps. To simulate many correlation times corresponding to asimulation of a few ms requires about 90 hours of IBM SP (SEABORG) time at 4×10-9
sec/particle-timestep. This is heroic, but the discharge time actually simulated is much smaller
than the plasma equilibration time.

The situation is similar for continuum codes. GYRO, in particular, has the most complete
physics to date for ion gyro-scale physics: ions and electrons (trapped and passing), magnetic
perturbations and collisions, real geometry, nearly full radius, finite ρ* = ρ ci a with profile andE×B shear stabilization and toroidal velocity-shear drive. Simulation of a single radial slice (∆r/a= 0.3) of DIII-D plasma for about 1ms takes five 24-hour submissions on 128 processors on
SEABORG. Scaling from DIII-D to ITER and taking full radius requires a factor of 15. The
1ms rescales to 2ms but it is still a long way from 3 sec confinement times. This reflects the factthat going beyond the previous state of the art from ions-only simulations to simulations withelectron dynamics requires an order of magnitude jump in computing power because therequired time step is at least 10 times smaller. Output from GYRO is shown in Fig. A8, showing
the strong anisotropy.

Fig. A8. Ion temperature-gradient turbulence in GYRO simulation. Extended structures along the field lines are evident.
[Courtesy J. Candy, GA]

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
18

D. NON-LOCAL PHENOMENA
We cannot leave this discussion without taking note of observed phenomena in tokamaks thatseem to defy the diffusion-convection picture. These include extremely rapid (i.e. faster than thetransport time scale) propagation of heat pulses (or “cold” pulses) resulting from sawtoothcrashes, impurity injection at the edge, or the H-mode transition itself. There are two current
approaches to explain these observations. First, there are empirical techniques based on theconcepts of self-organized criticality (SOC) and avalanches. Empirical models have given goodresults, but the underlying physics has not been elucidated. Second, the instabilities leading toturbulence usually have fairly well defined thresholds, e.g., a critical temperature gradient. A
rapid pulse may occur if the plasma is perturbed from its marginal condition. In this case, theunderlying physics is presumed known, but quantitative calculations have not been performed.
The simulation project must be prepared to incorporate developments in this area.
E. CHALLENGES FOR TURBULENCE SIMULATIONS
The bottom line is that we are currently far from the achievement of turbulence simulations withall the relevant physics on a scale suitable for integration with transport calculations. One
estimate, presented at the May 23 ISOFS Workshop, is that we are six orders of magnitude froma solution based on current computational methods and computers, and that we will make up
only four of these orders in the next few years by advances in computer technology and currently
envisioned schemes for exploiting the time-scale separation (see Sec. IXB). Among thecomputational and applied mathematical challenges are:••••Continuum kernels solve an advection/diffusion equation on a 5-D grid; we therefore
need: linear algebra and sparse matrix solvers (LAPAC, UMFPAC, BLAS), anddistributed array redistribution algorithms.
Particle-in-Cell kernels advance particles in a 5-D phase space and need: efficient“gather/scatter” algorithms that avoid cache conflicts and provide random access to field
quantities on a 3-D grid.
Continuum and Particle-in-Cell kernels perform elliptic solves on 3-D grids (often
mixing Fourier techniques with direct numerical solves).
Other Issues are portability between computational platforms, characterizing andimproving computational efficiency, distributed code development, and expanding theuser base.

Finally, we note that while core turbulent transport is extremely important (ability to predict
internal barriers, for example), edge turbulence, which has the same difficulties described here
and more (Sec. VIII), is critical, for the edge pedestal is the greatest source of uncertainty forreactor predictions.

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
19

IV. MACROSCOPIC SIMULATION
A. INTRODUCTION AND SUMMARY

Coupled magnetohydrodynamics (MHD) and Maxwell equations play a central role in modern
fusion plasma theory. First, MHD determines 2-D or 3-D magnetic equilibrium with nested
toroidal magnetic flux surfaces, which are crucial for magnetic plasma confinement. Next, MHD
describes both thresholds and nonlinear dynamics of device-scale plasma instabilities (so called
MHD instabilities). Very often these instabilities set the limits of the performance of fusion
devices.
Plasma dynamics can be completely described by the evolution of the distribution function
fα (r, v, t ) , for each particle species α, given by each species plasma kinetic equation, together
with the self-consistent evolution of the electric and magnetic fields, given by the Maxwell
equations, as given in Eqs. (A-1)-(A-6). Solving these equations for space scales and timescales
characteristic of large-scale instabilities in confined plasmas is computationally impractical. The
Extended MHD approach is to reduce the dimensionality of the problem, by multiplying the
kinetic equation by successive powers of the particle velocity v and integrating over velocity
space. If the underlying distribution functions have nice properties, such as a close-toMaxwellian velocity distribution, the resulting moment equations have fluid-like properties.
They are more tractable theoretically and computationally, although formidable problems may
still arise.
Magnetic fusion devices are very rich in MHD activity, some relatively benign, some leading tocatastrophic disruptions. Some of these are known as sawtooth oscillations, tearing and ballooning
instabilities, and resistive wall modes, whose general features and some quantitative predictions
can be given in terms of ideal or resistive MHD. As higher plasma temperatures are reached,more and more kinetic effects are required to be included.
B. EQUATIONS
On the time and space scales that the electrons and ions maintain local charge neutrality, thelowest order moment equations for the electron and ion species can be added together to form aset of equations for a plasma fluid with a density ρ = (Mi ni + me ne), fluid velocity v = (Mi ni vi +me ne ve) /( M i ni + me ne) ~ vi and pressure p = p i + p e. (Note that the symbol ρ is being usedhere for the mass density. The same symbol has elsewhere been used to represent charge density
and normalized minor radius.) The displacement current can be neglected in Ampere’s law (since∇ ⋅ J = 0), eliminating electromagnetic radiation and electrostatic oscillations. The MHD
equations can be written in a general form, in MKS units, as the low-frequency Maxwell
equations,∂B= −∇ × E ,∂t∇ × B = µ0 J ,∇⋅B=0 ,FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
20(A-13)

the continuity equation,∂ρ= −∇ ⋅ ρv ,∂t

(A-14)

 ∂v
ρ + v ⋅ ∇v = −∇ ⋅ P + J × B , ∂t(A-15)∂p5
2+ v ⋅ ∇p = − p∇ ⋅ v − (Π : ∇v − ∇ ⋅ q + Q).
33∂t(A-16)

the total momentum equation,and the energy equation,

The energy equation (A-16) assumes that the ratio of specific heats γ = γ e = γi = 5/3. In theseequations, µ0 is the permeability of free space, n is the number density, ρ is the mass density, v isthe center of mass velocity, B is the magnetic flux density, E is the electric field, J is the current
density, p is the scalar pressure, q is the heat flux, η is the electrical resistivity, the stress tensor is
P = pI + Π, where I is the unit tensor and Π is the traceless part of the stress tensor, and Q isother heat sources and sinks.
In general, the electron motion decouples from the bulk fluid motion, although the two are
related by v = ve + J/ne. For near-Maxwellian distribution functions, for example, decoupling
can occur due to the effects of a non-negligible ion Larmor radius (finite Larmor radius or FLR),
which is still small relative to the system size. The equation for the electric field, known as
Ohm’s law, comes from the electron momentum equation. Ignoring terms of order me/Mi, it is
E = −v × B + η J +1
1
J × B − ∇ ⋅ Pe −nene1  ∂v eρe+ ve • ∇v ene  ∂t(A-17)

It requires a pressure or temperature equation for the electrons,52∂pe+ ve ⋅ ∇pe = − pe ∇ ⋅ ve − (Πe :∇ve − ∇ ⋅ qe + Qe ) ,33∂t(A-18)

Additional detail, still within the confines of a two-fluid moment description, can be obtained bykeeping the anisotropies relative to the confining magnetic field, such as the two pressures p⊥ andp and/or the heat fluxes. The above equations then refer to the average quantities pj =(pj+2pj⊥)/3, etc.
To close the system, expressions for the higher order moments Π and q must be obtained
independently, from solutions to the kinetic equation. At high collisionality, these are the usualcollisional viscous stress tensor and the heat flux (proportional to the local velocity andtemperature gradients, respectively), and this leads to the Braginskii equations. At lower
FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
21

collisionality or long mean free path, these terms contain non-local kinetic effects. Proper closure
becomes a complex question that must take into account details of the confinement
configuration. (In toroidal systems, these non-local geometrical effects have been addressed in aflux-surface-averaged sense by neoclassical theory.) Unfortunately also, in this limit there is no
single, unambiguous way to define the set of “two-fluid” or FLR terms, so that models depend
upon a mixture of theoretical and practical considerations (see R. D. Hazeltine and J. D. Meiss,
Plasma confinement, Addison Wesley, 1992.)
Another approach that is being pursued is the so-called hybrid approach where a distribution ofparticles is used to provide closure to the fluid equations. There are 2 categories of the hybrid
approach that are being used: pressure coupling and current coupling. Even within thesecategories, there are several approaches.
In the pressure-coupling scheme, the distribution of particles is used to calculate the pressure
tensor, and then the velocity is advanced in time from Eq. (A-15). The pressure tensor has been
calculated differently by different researchers: from the gyroviscous stress tensor in terms ofgyrofluid moments; by using approximations to the gyroviscous stress tensor given in terms ofBessel functions in Fourier-Ballooning space; by using a particle Hall-MHD closure that uses
test particles to compute the off-diagonal elements of the pressure tensor; or by calculating thestress tensor directly from summing moments of Gyrokinetic particles.
In the current coupling scheme, the ion current is calculated directly from the particles, and no
fluid ion equations of motion is needed. This has been implemented using the full equation of
motion for the particles. This method is relatively straightforward but inefficient. There is also
discussion of implementing current coupling using gyrokinetic particles, in which case thepolarization current must be dealt with.
C. STATUS
Currently two state-of-the-art 3-D codes are being supported by the SciDAC Center for
Extended Magnetohydrodynamic Modeling (CEMM) project; M3D and NIMROD. These areboth focused on the modeling of linear and non-linear phenomena in fusion experiments that
require Extended-MHD descriptions.
The NIMROD code solves the primitive form of the plasma fluid-model in axisymmetrictoroidal, cylindrical, or periodic-linear geometry with arbitrary poloidal cross-sectional shape.
(The geometry must have an ignorable periodic coordinate, but the simulated dynamics are fully
three-dimensional.) The user selects which terms are retained in Ohm’s law through an input
parameter. The semi-implicit numerical method is used to advance the solution from initialconditions. This avoids severe time step restrictions associated with wave-like normal modes ofthe system, sound, Alfvén, and whistler waves—while avoiding numerical dissipation. For
accuracy at time steps that are orders of magnitude larger than explicit stability limits, the semiimplicit operator for mass motions is based on the linearized ideal MHD energy integral. Matrix
inversion is accomplished by parallel preconditioned Krylov methods, which is the mostcomputationally demanding part of the time advance. Performance is therefore dependent on theeffectiveness of the preconditioner.

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
22

The spatial representation of NIMROD is an important feature of the code. NIMROD uses a
combination of logically quadrilateral and triangular finite elements for the poloidal plane andpseudospectral collocation for the periodic direction. The polynomials used for finite element
basis functions are selected by the user for optimal efficiency, and poloidal mesh lines need not
be orthogonal. For many fusion problems, accuracy is improved by aligning grid lines with theequilibrium flux surfaces inside the separatrix. The grid can also be packed around low orderrational flux surfaces to efficiently resolve the small spatial scales that arise at high Lundquist
number S. Triangular meshing outside the last closed flux surface allows complicated, realistic
boundary shapes.
The M3D code is a parallel code that is especially suited for geometries with inherently threedimensional boundaries, e.g. stellarators, but can also be used to simulate axisymmetric devices.
M3D consists of two parts, a mesh module and a physics module. The mesh module contains
the grid, implementation of differential and integral operators, I/O, and inter-processor
communication. The physics module handles time advancement of the equations and contains ahierarchy of physics levels that can be invoked to resolve increasingly complete phase-spaces, andtherefore provide increasing realism. The module includes resistive MHD, two-fluid, and kinetic
particles. Electrons are represented as a fluid with an approximate fluid closure. M3D uses a
stream function/potential representation for the magnetic vector potential and velocity that has
been designed to minimize spectral pollution. Parallel thermal conduction is simulated with the“artificial sound” method. The solution algorithm is quasi-implicit in that only the most timestep limiting terms including the compressional Alfvén wave and field diffusion terms are
implemented implicitly, with explicit time stepping used for the remaining terms.
A three-dimensional mesh is utilized to facilitate the resolution of multi-scale spatial structures,such as reconnection layers and to accommodate fully three-dimensional boundary conditions
that occur in stellarators or the evolving free boundary of a tokamak bounded by a separatrix.
The mesh uses unstructured, 3-D piecewise-linear triangular finite elements in the poloidal
sections. The domain decomposition consists of slicing the toroidal geometry into a set ofpoloidal planes with each poloidal plane further partitioned into equal area patches. One or moreof the poloidal patches are assigned to each processor. The fluid part of each time step consists of
uncoupled 2-D scalar elliptic equations that are solved concurrently within each poloidal plane.
The PETSc library has been used extensively to provide a portable, efficient parallel
implementation for the elliptic equations that need solution at each time step. These are solved
with a Krylov accelerated iterative scheme that uses the overlapped Schwarz method forpreconditioning. This leads to excellent parallel scalability.
D. CHALLENGES
These codes require high resolution and many time steps to give an accurate representation ofmodern fusion experiments. To see this, recall that resistivity effects are characterized by aresistive diffusion time scale, τ R ~ µ0L2 /η , which is much larger than the Alfvén-wave transit
time, τ A ~ L ρµ 0 / B . (Here L ~ 1 m is the spatial scale of the device). In fusion machines the
Lundquist number, S = τ R /τ A , is of the order of 10 8 . As a result, even though resistive effects
determine the physics of the process, they actually become important only within some very

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
23

narrow layers. Proper resolution of these layers as well as strong anisotropy of plasma properties
along and across the magnetic field imposes serious computational challenges.
We can estimate the computational resources required to carry out the necessary simulations forparameters typical of present and proposed experiments as shown in Table A-I.
TABLE A-I. Typical dimensionless parameters for present and proposed experiments
ITER
parameter
name
CDXU NSTX CMOD DIII-D FIRER (m)Te [keV]βS1/2(ρ*) –1 =Ba/T1/2
a/λe

Major radius
Elec. Temp.
Plasma beta
Inv resis. length
Ion number
Recip. norm
elec skin depth

0.3
0.1
0.01
200
40
250

0.8
1.0
0.15
2600
60
500

0.6
2.0
0.02
3000
400
1000

1.6
2.0
0.04
6000
250
1000

2.0
10
0.02
20,000
500
1500

5.0
10
0.02
60,0001200
3000

Estimates based on explicit time-stepping with no grid refinement. Let us first estimate thecomputational requirements for a 3-D calculation with uniform zoning of size ∆x and an explicit
time-stepping scheme based on the CFL criteria for the poloidal Alfvén wave, i.e. ∆t = ∆x / VAP.
For a 3-D mesh of linear dimension N, i.e., N 3 mesh total mesh points, it would take N time
steps to calculate one Alfvén wave transit time τA = a / V AP. Typical ideal and resistive MHD
instabilities would grow on the timescales TIDEAL ~ β-1/2 τ A and T RESIS ~ S1/2 τ A , requiring about
β-1/2 N and S1/2 N time-steps, respectively. Thus, the total number of space-time points required
to compute an ideal or resistive instability would be about β-1/2 N4 (ideal) and S1/2 N4 (resistive).
The plasma beta, denoted by β, is the ratio of plasma pressure to magnetic pressure. Its value isan important measure of confinement quality.
Current experience shows that with real performance of about 100 Mflops/processor and oforder 1000 processor-hours, we can compute a problem with 1003 mesh points for 104 time
steps, using a complex fluid model with the compressional wave and field terms implemented
implicitly. This is about 1010 space-time points in 3 × 1014 operations. This is easily sufficient
resolution and time steps to calculate an ideal mode in CDX-U, and is nearly sufficient to study
the initial growth phase of a nonlinear tearing mode. We see this since the number of linear
mesh points is comparable to: the linear tearing-layer width, S −1 2 , the ratio of the system size tothe ion Larmor radius 1 ρ *, and the ratio of the system and the ratio of the system size to theelectron collisionless skin depth, a λe . These are the relevant lengths that enter the two-fluid
Extended MHD model.
The question is: what type of computer power is needed to study this physics in a larger, hotter
device with a stronger magnetic field? We see from Table A-I that depending on what scale
length needs to be resolved, the number of mesh points in a linear direction will increase by
about an order of magnitude as we go from CDX-U to DIII-D. The increase in the total numberof space-time points would be the fourth power of this factor times another scaling factor thatwill between unity (for ideal scaling) to about 10 (for resistive scaling). Thus, the number of
space-time points required would increase anywhere from 104 to 105. Running on a 10 Teraflops
FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
24

(delivered) computer for 3 days would correspond to about 3 × 1018 floating point operations
which would be about 104 times greater than what was quoted above for what is available to us
today, so a full DIII-D calculation might be feasible with this hardware increase alone.
Grid refinement, implicit time stepping and improved algorithms. It is straightforward to see thatthe above scaling estimates can be gross overestimates if we take into account improved
algorithms and meshing. For example, for a field-line following mesh and with adaptive meshrefinement, the total number of mesh points should only have to grow linearly as we go to larger
machines and higher resolution, rather than cubic. With implicit time differencing, the time step
will not have to decrease nearly as fast as linear with zone size. More efficient solvers can give anadditional factor. Thus, with these computational improvements, some of which have alreadybeen implemented to varying degree in M3D and NIMROD, we can realistically expect to beable to calculate modes using Extended MHD models in DIII-D, NSTX, and CMOD in thetime frame of this project, and even FIRE and ITER calculations might be within reach. Of
course, it also may not be necessary to simulate the exact parameters of a machine if we candetermine scaling relations from doing a series of calculations at reduced parameters.

V. PLASMA EDGE PHYSICS AND PLASMA-WALL INTERACTION
A. INTRODUCTION

The edge plasma, which bridges the hot plasma core and the material wall (see Fig. A9) plays acrucial role in both overall plasma confinement and plasma-wall interactions. Some examples of
the impact of the edge plasma are:
1. Changes in edge plasma parameters can lead to dramatic improvement in core plasma
confinement in the H-mode via the formation of transport barriers, regions with reduced
plasma turbulence. The height of the pedestal of this barrier plays a crucial role in theperformance of a fusion reactor;
2. Practically all magnetic fusion devices have a density limit, which may be due to strong
anomalous plasma transport at the edge;
3. There is a strong, but poorly understood, influence of neutrals and wall conditions on
plasma confinement, e.g., one of the best shots in the TFTR tokamak was enabled bylithium conditioning of the walls;
4. Heat load on first wall, which is determined by edge plasma conditions, is a serious issue
for reactor-relevant conditions;
5. Wall sputtering, transport of ions and neutrals, including hydrocarbons, in the edge
plasma, and deposition processes caused significant accumulation of tritium in the firstwall of both TFTR and JET tokamaks.

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
25

separatrix

SOL

core edge
first wallFig. A9. Schematic view of edge plasma region in tokamak

The plasma edge region has many of the same problems as the core, but it also has a number of
attributes that make it crucially distinct from the core. In particular, in toroidal magnetic fusion
devices, such as tokamaks, field lines in the core of the device lie, at least approximately, onclosed toroidally shaped flux surfaces, giving rise to good confinement of particles and heat.
However, on the periphery of such devices, there inevitably exists a scrape-off layer (SOL), where
magnetic field lines are in direct contact with material surfaces. Due to competition of paralleland perpendicular plasma transport, two-dimensional effects are strong in the SOL. In addition,so-called divertor configurations are common to most of the large toroidal devices, in which extra
field coils are added to make the magnetic field intersect material surfaces at a location relatively
remote from the main plasma volume, as shown in Fig. A9. In divertor configurations themagnetic separatrix divides closed magnetic field lines in the core-edge region from open ones in
the SOL.
In the core plasma, the separation of scales is both a blessing and a curse. On the one hand, as ingyrokinetic theory, an ordering exists that allows one to solve the equations by averaging oversmaller scales and ignoring variations in the longer scales. On the other hand, the same scale
separation contributes to the difficulty of a first-principles whole-device simulation. In the edge,
by contrast, the spatial scales are compressed, tending to invalidate the ordering schemes.
Gyrokinetic ordering, for example, can break down, i.e. ρ* ~ 1, and the collisionality varies along
the field from the long mean-free-path to the short mean-free-path regime.
Thus, the solution of the fundamental equations becomes much more difficult, but if solutions
are obtained, they should be readily assimilated into the global simulation–or become amicrocosm of the global simulation. Indeed, one of the focused integration initiatives (FIIs)
might be devoted to this entire region.
B. PROPERTIES OF THE EDGE PLASMA
Plasma near the separatrix, being impacted by fast parallel transport to material surfaces, tends tohave steep radial gradients in temperature and density and to be relatively cold. The lowtemperature, coupled with proximity to bounding surfaces, results in a relatively highconcentration of neutral gas and impurities. These properties lead to a relatively strong role foratomic physics processes: ionization, recombination, excitation, and radiative transport.
FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
26

There is general consensus, supported by both analytic theory and numerical simulations ofanomalous transport, that the poloidal E×B flows, spontaneously generated by the nonlinear
dynamics, play a central role in regulating the saturation level of the turbulence and the resulting
cross-field transport. Strong shear of this poloidal flow tears apart convective eddies and reduces
the turbulence level and cross-field transport by forming the H-mode transport barrier. Due tothe interplay between significantly different physics on open and closed magnetic flux surfaces,
strong shear of the radial electric field and E×B poloidal flow arises, and, as a result, an H-mode
transport barrier may be formed somewhat inside the separatrix. Reduction of anomaloustransport at the barrier causes steepening of the plasma temperature and density profiles in thisregion. As a result, strong and repetitive MHD modes can develop causing ELM bursts. It isbelieved that the MHD modes responsible for ELMs are the so-called ballooning and peelingmodes driven by plasma pressure and electric current gradients.
With increasing plasma density, the plasma particle flux to the divertor targets starts to decreaseas the detached divertor regime is being formed. It is rather well understood and shown with both
simplified analytic models and sophisticated 2-D plasma transport codes such as UEDGE thatplasma-neutral coupling and atomic-physics effects, including impurity radiation and plasma
recombination, play key roles in establishing this regime. However, there are strong indications
that cross-field plasma transport also plays an important role here.
Disruption of the discharge for densities above the density limit is believed to be due to the mixed
effects of enhanced plasma transport to the first wall and thermal collapse of the plasma due toimpurity and hydrogen radiation. Interestingly, both the DBM (Univ. Maryland) and BOUT
(Lawrence Livermore National Lab.) edge-plasma turbulence codes show trends in plasmatransport enhancement at high plasma densities that are somewhat similar to that seen inexperiments.
We note that in both improved confinement H-mode and standard L-mode there is rather
strong interaction of plasma with first wall material surfaces resulting in both sputtering and (reco-) deposition. These are rather complex processes involving ion transport and neutral-particle
transport in edge plasmas, chemistry associated with both heavy particle interactions andinteractions with electrons, and surface effects (implantation, collision cascades, deposition,
adsorption, desorption, diffusion, etc.). In all regimes the first wall is a huge reservoir ofhydrogen isotopes, which plays a dominant role in neutral gas recycling.
C. COMPUTATIONAL CHALLENGES AND CODES FOR THE EDGE REGION
Fusion plasmas in general have a large span of spatial and temporal scales, from the electron
gyroradius and cyclotron frequency to the device size and the energy confinement time. Thislarge span makes simulation a substantial challenge. In the hot core plasma, there is often a wide
separation between the space and time scales characterizing turbulent fluctuations and thosecharacterizing evolution of the equilibrium. This scale separation facilitates the use of separate
simulation codes to describe turbulence and transport. However, in the edge region, this scale
separation can become small or even non-existent, leading to the challenge of combining a wide
range of turbulence and transport scales into a single, large-scale simulation. In addition to theplasma-physics scales, the presence of important atomic physics processes introduces new length
FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
27

and time scales, which create a range larger than that for the core. Ionization, recombination,and charge-exchange rate coefficients for hydrogen and impurities are in the range of 10-14 - 10-13
m3/s, which for typical densities of 1020 m-3 give time scales of 10-6 - 10-5 s. On the other hand,near-unity recycling from the hydrogen-saturated material surfaces, where each incident ion
results in a neutral hydrogen atom injected back into the edge plasma, yields a long time scale,
10-1 s, for establishment of equilibrium profiles.
The edge region of fusion devices generally has substantially lower plasma temperatures than thecore, resulting in Coulomb collisional mean-free-paths parallel to B being much less than theconnection length, i.e., the parallel distance traveled in making one poloidal revolution.
Furthermore, the gradient scale-lengths and observed turbulent wavelengths perpendicular to B
are usually greater than the ion gyroradius. Because of the spatial localization provided bycollisions and B, fluid models have been adopted to give the basic description for both
turbulence and transport simulations. However, there is growing concern in the community thatfluid models are not adequate for plasma conditions in the H-mode pedestal region, andtherefore they should be replaced with more accurate, but much more complex, kinetic models.
The principal edge-plasma turbulence codes within the U.S. community are presently BOUT
and DBM. These are 3-D turbulence codes dealing with a fluid plasma description based on theelectromagnetic Braginskii equations for plasma vorticity, density, electron and ion temperaturesand parallel momentum. The BOUT code is non-local, can describe a magnetic X-point
geometry on both sides of the separatrix, and is based on a toroidal segment simulation volume,
while the DBM code is a flux-tube code (presently without an X-point) and well suited toparametric studies. With sources added in the core-edge region and sinks in the SOL, BOUT has
begun to follow some short-time profile evolution in response to the turbulence.
The numerical algorithms used in BOUT consist of finite-difference equations in 3-D where theresulting ordinary differential equations for the time dependence of each cell variable areadvanced with the fully implicit Newton-Krylov solver PVODE. BOUT is written in C. The
implicit integration increases the time step by a factor of 3-6 without preconditioning.
Parallelization is obtained by domain decomposition in the direction along B, utilizing MPI for
message passing. Because of weak coupling between the domains, scaling with the number ofprocessors is essentially linear up to 120, but using 64 processors is more typical, and simulations
have been done on various parallel platforms (IBM SP, T3E, SUN, DEC and Linux clusters).
The DBM code is written in Fortran 90 and solves the reduced Braginskii equations in a general
flux-tube magnetic geometry, without (as yet) a magnetic X-point. It utilizes fourth-order spatial
finite differencing and a second-order-accurate trapezoidal leapfrog time advance. The
communication routines are MPI-based and at present, like BOUT, involve 1-D domain
decomposition.
In addition, the edge plasma community has extensive experience with kinetic neutral transportand related atomic physics via DEGAS-2 (PPPL) and TNG (UCSD), as well as the 2-D coupled
fluid plasma/neutral transport via UEDGE (LLNL), all of which provide an excellent base forplanned extensions. UEDGE is an implicit time-dependent code capable of very long-time
simulations of profile evolution that includes the important neutral particle sources fromrecycling surfaces and gas puffing, plus impurity species.
FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
28

However, even though there is a consensus that neutrals can be an important ingredient in bothedge-plasma turbulence and in the formation of the H-mode barrier, so far edge plasma
turbulence codes are lacking the proper physics that would take these effects into consideration.
To describe erosion of the first wall under normal operation conditions the REDEP/WBC
package is often used in the U.S. edge plasma community. The WBC is a Monte Carlo code,
which computes impurity-atom and ion motion at the kinetic description level, including
sputtering and deposition processes, and both elastic and inelastic impurity-plasma collisions.
The code is very time consuming to run; therefore it is often run separately, and its output used
as input to the REDEP code, which uses cruder models for impurity transport. As an input forplasma parameters the REDEP/WBC package uses either experimental data or UEDGE
modeling results. We note that the cross-field impurity-ion transport built into this package israther rudimentary

VI. EQUILIBRIUM EVOLUTION AND CONFINEMENT
A. INTRODUCTION

P lasma transport or confinement will be a major component of the integrated simulation
project. Processes described in the preceding section determine the sources and sinks that drive
the plasma fueling, heating, and rotation. In the following section we discuss the confinement ortransport mechanisms that balance these. First, however, let us set up the framework in whichthese operate. With exceptions to be discussed below, the plasma equilibrium consists of nested
flux surfaces. Because plasma parameters equilibrate rapidly along the magnetic field, they arenearly constant on these surfaces, and thus transport takes place principally across magnetic
surfaces, and can be described in one spatial dimension, denoted by a generalized radius or flux
label ρ. The term 1-1/2 D transport refers to the fact that ρ must be related geometrically to the2-D equilibrium. These 1-1/2 D transport codes are central to today’s efforts at integrated
simulation. We will come back to this point later.

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
29

Fig. A10 PIES equilibrium of a low aspect-ratio stellarator [Courtesy J. Lyon and W. Houlberg, ORNL].

The breakdown of the above argument gives rise to a physics and computational challenge thatwill be central to the integration initiative. Islands and ergodic regions can form, owing to MHD
instabilities and field errors in tokamaks, as well as in 3-D equilibria in stellarators. Figure A10
depicts a cross-section of a low aspect-ratio stellarator equilibrium showing islands and ergodic
regions. (This is a puncture plot with each dot representing the passage of a magnetic field line.)
These regions evolve on the transport time scale.
B. TRANSPORT EQUATIONS
Despite these caveats, tokamak transport is reliably described by 1-D surface-averaged transport
equations in most circumstances. For stellarators, a design goal is to minimize the externallydriven islands and stochastic regions so that similar transport studies are approximately valid.
Construction of a set of 1-D equations makes use of characteristic time scale separation. The
fastest or Alfvén time scale (~µs) is assumed to establish the basic magnetic geometry. This iscomputed with MHD equilibrium codes that are free from the Alfvén time scale so this timescaleis effectively eliminated from the problem. Likewise, particle densities and temperatures are
assumed to equilibrate rapidly along the magnetic field lines, so that they can be considered 1-D
functions of the flux coordinate function, thus eliminating the fast parallel transport time scalefrom the problem. The 1-D functions evolve on timescales characteristic of cross-field transport,
which is one of the slowest timescales under consideration. This timescale establishes the profiles

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
30

of thermodynamic quantities, density, temperature, and angular momentum. Careful selectionof variables is required to take advantage of this timescale separation.
The geometry of the flux surfaces is specified by solving the MHD equilibrium equations:
J × B = ∇p , ∇ × B = µ0 J , and ∇ ⋅ B = 0 . In two-dimensional geometry, appropriate fortokamaks without magnetic islands, flux surfaces are known to exist, and this leads to a nonlinear
partial differential equation, the Grad-Shafranov equation, for the poloidal magnetic flux Ψ as afunction of the cylindrical coordinates R, Z in the poloidal plane. Different techniques are used
for fixed boundary calculations where the plasma/boundary interface is specified, and freeboundary calculations where the plasma/boundary interface is determined self-consistently fromthe actual coils that produce the confining magnetic field. In three dimensions, there is anadditional complexity in that there is no guarantee that the nested flux surfaces exist. The
VMEC 3-D equilibrium code assumes the existence of these surfaces, whereas the PIES 3-D
equilibrium code does not.
Figure A5 (or A10), without the islands or ergodic regions, is a typical tokamak (or stellarator)
result. The radial variable ρ(Ψ) is usually defined in terms of the plasma volume or magnetic flux
(integral of B ⋅ eˆ over a surface perpendicular to eˆ ). For example, ρ = a( Ψ Ψa ) ,
1/ 2

1/ 2

ρ = a(V /Vtot ) or ρ = a(Φ /Φtot ) , where a is the nominal minor radius of the device, V is thevolume and Φ is toroidal flux within a flux surface.

The transport equations themselves are obtained by taking velocity moments of the Boltzmann
equation (A-1) and computing appropriate flux-surface averages,

A =

2π
2π
1
dφ ∫ dθ g( ρ,θ, φ ) A( ρ,θ, φ ) ,∫00V ′( ρ )2π

2π

V ′( ρ) = ∫ 0 dφ ∫ 0 dθ g( ρ,θ, φ ) ,(A-19)(A-20)

where V ′ = dV dρ , g is the Jacobian of the coordinate transformation and A is the function to
be averaged. The resulting time-dependent continuity and energy transport equations are of the
form:
1 ∂
∂n a
=−
(V ′ Γa ⋅ ∇ρ ) + S pa ,∂t
V ′ ∂ρ
3 ∂( n a Ta )1 3
∂  Γ ⋅ ∇ρ 
= − V ′ Q a + Ta Γa  ⋅ ∇ρ  − n a Ta 


2 ∂t
2
∂ρ  n a 
V ′ − (∇ua ) : Π a + ∑ ∆Qaa ′ + SEa ,(A-21)(A-22)a′whereΓa ⋅ ∇ρ and Q a ⋅ ∇ρ are transport fluxes, S pa and SEa are sources of particles andenergy, (∇ua ) : Π a is the viscous heating, and Qaa ′ is the energy exchange between species.
There are similar equations for momentum balance. Additional constraints are that fast energy

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
31

exchange between ions forces their temperatures to be nearly equal and that quasi-neutrality,∇ ⋅ J = 0, determines ambipolarity constraints,

n e = ∑ Z a n a , Γe = ∑ Z a Γaa= ions(A-23)

a= ions

The fluxes, viscous heating, and sources are inputs to the code and are discussed in more detail inthe preceding and following Appendices. Kinetic theory is used to provide closure, whereby
Γa , Q a , and Π a are expressed in terms of the thermodynamic variables n a , Ta and their
gradients. For example, the particle flux is usually expressed in the diffusive-convective form∂n a
2 1/ 2+ ∇ρ
n a una ,∂ρ∂T2
2 1/ 2Q a ⋅ ∇ρ = − ∇ρ n a χ a a + ∇ρn a Ta uqa∂ρΓa ⋅ ∇ρ = − ∇ρ

where ∇ρ2

2

Da

(A-24)(A-25)

is a geometrical factor, Da and χ a are the particle and heat diffusivities, and una

and uqa are the convective velocities. The latter may themselves be represented in terms of otherplasma gradients and referred to as off-diagonal elements.
The magnetic fluxes also evolve. This is typically described in a Grad-Hogan scheme, in which
the toroidal flux is taken as the reference frame and is updated with infrequent calls to the MHD
equilibrium equation. The poloidal flux evolves relative to the toroidal flux and is described byresistive diffusion, making use of Faraday’s law (see Eq. (A-3)) and Ohm’s law in a flux function
form that relates the toroidal current density as the secondary of a transformer to the externalloop voltage. The time scale is governed by the parallel electrical resistivity.
C. USES OF TRANSPORT CODES
Codes that solve the transport equations as described here form the core of current state-of-theart integrated modeling. These codes have proliferated over the years to serve a multitude of
purposes, and each achieves good results within its range of validity. We distinguish two basic
modes of operation. In the interpretive mode, all (available) plasma-parameter profiles and theirtime derivatives are inferred from experimental data, including the equilibrium geometry, whichis calculated in an MHD code using magnetic data as constraints. The sources, such as energy
deposition, are likewise calculated using the measured profiles. The code then solves (A-21) and(A-22) to calculate the transport fluxes or diffusivities, Eqs. (A-24)-(A-25), which are used todetermine empirical confinement scalings and the like. In the predictive mode, on the other hand,models for the diffusivities or fluxes are employed, and the code calculates the time evolution of
Eqs. (A-21) (A-22) to determine the plasma profiles, which can then be compared withexperiment.
The major interpretive codes in use today, TRANSP and ONETWO, contain detailed modelsof all the principal sources and sinks that are relevant to tokamaks. Each can be run,alternatively, in the predictive mode. Indeed, in the absence of one or more pieces of

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
32

information, a subset, such as the current or poloidal-field evolution, may be simulated within
an interpretive run.
Other codes, such as BALDUR, WHIST, CORSICA, or TSC, are simulation codes whose mainuse is testing transport models or performing design studies. These carry varying levels of detailoutside the particular feature being tested. For example, a core thermal-transport model may betested by running only the temperature evolution equation, holding the density fixed, setting aboundary condition just inside the edge, and taking the source as calculated in an experimental
interpretive run. The goal would be to test the temperature profile evolution withoutcomplicating distractions. Similar simulations in which all of the density, temperature, andmomentum equations for each species plus current penetration could be run, while still using
simple boundary conditions. At the other end of the spectrum, one might wish to explore reactor
designs, studying wall conditions, connection to external circuits etc., while employing a very
simple plasma model based on empirical scaling. Each of these problems requires different
functionality in the code to make efficient use of the available computer time.
D. CHALLENGES OF TRANSPORT MODELING
Flexibility. A principal challenge of the simulation initiative is to preserve present functionality
within a flexible configuration that allows the user to select from the full range of options andefficiently run his/her case, from the very simplest model to the most detailed full integration.
Rotation and radial electric field. The radial electric field, poloidal and toroidal rotation, and thepressure gradient are all coupled through the radial force balance. Both the physics and thecomputation of these processes require continued development and will form an essential aspectof integrated simulation. For example, bifurcations are observed experimentally, in which theplasma jumps suddenly to an alternate state, often forming a transport barrier. Enhanced
diagnostics and developing theoretical models imply scenarios in which sheared E × B flow
suppresses turbulence, while turbulence can drive locally sheared poloidal rotation. Also,∂E ρ ∂ρ squeezes orbits and reduces collisional transport, while toroidal rotation enhances
MHD stability. Toroidal rotation generally appears diffusive and governed by turbulent
transport but can appear even in the absence of apparent torque. Other issues are whether weshould treat E ρ or poloidal rotation as the independent variable, and whether we can express
bifurcation threshold conditions in terms of macroscopic quantities. Internal barriers are strongly
localized and move in time, indicating the need for dynamic gridding algorithms for better
resolution. Finally, in stellarators, the ambipolarity constraints are expected to govern E ρ , and intokamaks, broken toroidal symmetry can do the same.
Transport in the edge and SOL. The basic equations in the edge and SOL are much the same as in
the core, but the computational challenges are even greater. Physics and computational
challenges include formulation of both turbulent and collisional transport models in steep
gradient regions, formulation of the coupling/transition between the closed field lines (slowradial transport timescale) and open field lines (faster parallel timescale), and accommodation ofmultiple timescale phenomena e.g., edge localized modes (ELMs) of the H-mode layer. (See Sec.
VIII.)FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
33

Non-nested surfaces. Finally, one of the most difficult challenges of integrated core modeling is toextend all the considerations of this and the adjacent two sections to take into account theislands and stochastic regions illustrated in Fig. A10 .

VII. EXTERNAL SOURCES
A. INTRODUCTION

E xternal systems that add mass, momentum, or energy to a plasma are essential tools in the
successful efforts to obtain good performance both in present experiments and in any future
experiments and reactors. At present, the external sources include beams of neutral atoms thatcan carry energy, particles, and angular momentum across magnetic fields; radio frequency waves(RF) whose interactions with a plasma can be used for heating, current drive or flow drive; highspeed pellets of frozen fuel gas that can deliver particles deep into the plasma core; and gas
fueling which supplies particles to the plasma edge. In a burning plasma experiment or reactor,
the fusion-produced alpha particles, which have some properties in common with beam-injectedparticles, will provide the principal energy source.
External sources are important, not only to provide fuel (D and T in a reactor) and energy(heating the plasma to 10 keV or more), but also to provide essential elements of control. The
plasma current, pressure, and rotation-velocity profiles all are crucial to plasma performance andcan be influenced by these sources. We also have sources such as impurities and gas from plasmafacing components, which can have undesirable consequences, as well as losses from radiation.
The computational capabilities related to neutral beam injection are very well developed and arepresently integrated into many codes. The work could be easily ported to codes developed in
new initiatives. Wave-plasma interactions at RF frequencies (from ion to electron cyclotron
frequencies) are the subject of intense ongoing research including fusion SciDAC activity.
However, the scope of needed work, and the ability to provide interactive coupling with otherplasma codes, extend far beyond the SciDAC activity. This area includes many distinct problems
and will be an essential element of future plasma prediction, interpretation, and control schemes.
Work on the ablation and subsequent transport and deposition of fuel from injected pellets is ina relatively early phase of development but certainly amenable to computation. However thephysics of fueling in general is not in a satisfactory state at present and would benefit from basic
theory studies of particle transport.
B. NEUTRAL BEAM INJECTION
A number of well understood physics processes are involved in neutral beam injection (NBI) at
energies in the neighborhood of 100 keV. The capture of neutral atoms by ionization and charge
exchange is modeled by Monte Carlo calculations. A thermal neutral atom formed by charge
exchange may escape or be re-ionized. The orbits of fast charged particles are followed and theirslowing down and scattering are computed in a Fokker-Planck model solved by Monte Carlo.
Secondary processes such as re-neutralization of energetic ions by charge exchange are also
calculated.
FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
34

For tokamaks the TRANSP neutral injection package is the most widely used. It has recently
been extracted from the integrated interpretive modeling code TRANSP as a separate module for
the National Transport Code Collaboration (NTCC). See the web site
http://w3.pppl.gov/NTCC/NUBEAM/UserGuide. Much of interface code is automaticallygenerated by Python script. Input and output is passed in f90-derived data types, employing
internal 2-D flux coordinates (ρ , θ ), where ρ is the normalized minor radius, and θ is thepoloidal angle (the short way around in Fig. III.1 of the main report). Outlines for conversion
among other common coordinate representations are included. Inputs consist of MHD
equilibrium quantities (magnetic field geometry, plasma pressure in space), profiles of density
and temperature of a large number of particle species, descriptions of injected neutral beamgeometry and characteristics, atomic physics data and cross sections, wall geometry, and controls
for code operation and inputs to auxiliary models. Outputs are radial profiles of power
deposition, driven current, ion/electron source rate, and rotation, 2-D profiles of neutral and fastion density, and the fast ion distribution in four phase-space dimensions (ρ, θ, E, v||/v).
The issues for this kind of modeling are more computational than physics: increasing speed,
increasing resolution in space and velocity space, and extension to 3-D equilibria for stellaratorsand other non-axisymmetric configurations.
An entire separate field of study, however, concerns possible instabilities driven by populations ofhigh-energy particles. These are seen experimentally and considerable progress has been made in
explaining the observations by theory and computation. These are important because they canlead to anomalous losses of fast particles and reduce heating efficiency.
Alpha particles produced at 3.5 MeV in D-T fusion reactions have many of the properties ofhigh-energy beam particles, resulting in alpha heating and including the prospect of driving
instabilities, and can be modeled similarly. The chief differences are the much higher alphaenergies and their isotropic initial distribution. Because large numbers of alphas will not be
produced until a burning plasma is actually achieved, it is extremely important to include thisphysics in our integrated computation and in this way anticipate burning plasma performance.
C. EDGE PARTICLE SOURCES
Modeling of edge particle sources is a discipline unto itself and must be included in the edge
plasma studies. (See Sec VIII.) Edge conditions in fusion devices have a tremendous influence on
the bulk plasma behavior. Neutral particles in the edge also influence other physics processes
besides fueling. For example, charge exchange of plasma ions with edge neutrals exerts torqueaffecting plasma rotation, and neutrals participate in edge instabilities. Edge particle sources are
inherently three-dimensional. Several modeling codes presently exist. Most are completely 3-D
Monte Carlo (EIRENE, DEGAS). A fast neutral code (NUT), based on zonal integration has
recently been placed in the NTCC module library.

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
35

D. RADIO-FREQUENCY WAVE HEATING AND CURRENT DRIVE
Besides inductive Ohmic heating, common to all tokamaks (but usually absent from stellarators),
there are three important wave frequency regimes to consider, each associated with a particle orplasma resonance. These are:• Electron cyclotron range of frequencies (ECRF), ω ≈ ω ce .
f ~ 100 GHz ( τ ~ 10-11 sec), λ ~ 0.3 cm, where λ is the wavelength. The launcher is far
from the plasma; the waves propagate in free space and can be computed withgeometrical optics.
• Lower hybrid range of frequencies (LH), ω ce << ω << ω ci .
f ~ 0.5 - 5 GHz ( τ ~ 10-10 sec), λ ~ 1 cm. The launcher is near the plasma; the waves do
not propagate in vacuum but are usually computed with geometrical optics.
• Ion cyclotron range of frequencies (ICRF), ω ≈ ω ci
f ~ 100 MHz ( τ ~ 10-8 sec), λ ~ 10 cm. The launcher is near the plasma, the waves do
not propagate in vacuum, and they require solution of wave equation. Alfvén waves atfrequencies ω << ω ci have also been considered for heating and current drive.
Waves can have strong interactions at localized regions in space due to various kinds ofresonance. Either the wave velocity matches the particle velocity v wave = ω k = v particle , and theparticle sees a steady, accelerating electric field, or the wave frequency in a frame moving withthe particle (Doppler effect) matches a harmonic of the particle cyclotron frequency,
ω − kv particle = lω c , and the gyrating particle sees a steady component of the electric field. The
effect is an energy and momentum kick each time the particle passes through resonance. These
kicks accumulate over time to produce energy gain and directed velocity resulting in bulk plasma
heating, very energetic tail populations, electric current and fluid drift or flow. In Fig. A5, wehave shown a tokamak cross-section with schematic launchers, a wave, and a resonance layer.
(Also depicted is a chain of magnetic islands. See Sec. IV.) Similar resonances also occur
throughout the plasma when spontaneous turbulence arises.
This illustrates an important difficulty in solving the basic equations, Eqs. (A-1)-(A-3), namelythat for perturbations from equilibrium, the distribution function f a ( x,v,t ) can be highly
structured, requiring very fine resolution in some regions of both configuration and velocity space.
Let us consider the full-wave calculations as they apply, for example to ICRF heating at
frequency ω. The wave propagation and plasma response are described by Eqs. (A-1)-(A-5) withdisplacement current included. We have now returned full circle to the highest frequency rangesof the whole integrated simulation. For this problem, we exploit the separation of time scales andthe small amplitude of the waves to linearize the fields B(r,t ) = B 0 (r ) + B1 (r )e−iωt ,E(r,t ) = E1 (r )e−iωt , and apply the quasi-linear method of solution to the plasma response foreach species, f a = f 0 + f1e−iωt + f 2 , where f1 is the fast response producing the wave current J in
Eq. (4.4), and f 2 is the slow time scale response that gives the power deposition and equilibrium
evolution. We give a schematic representation of the equations, in which the fast part of Eq. (A1) reduces to

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
36

iωf1 + L1{ f1} = L 2 {E1, f 0 } + L 3 {B1, f 0 } ,(A-26)

where L1-L3 are complicated linear operators. The slow background response then is the solutionto∂f 2+ L1{ f 2 } = Q1{ f1,E1} + Q 2 { f1,B1} + C( f 2 )∂t,(A-27)
where Q1 and Q2 are quadratic or bi-linear operators, proportional to f1E1 and f1B1 , where only
the “zero-frequency” components are retained. The set is completed with boundary conditions
that connect to the vacuum region and antennas that drive currents J ant e−iωt . To solve theseequations, a spectral representation is employed in which the plasma is treated as locally uniformwith Fourier expansions in three dimensions.
The wave code with the least restrictive approximations is the All Orders Spectral Algorithm(AORSA). Discretizing using the method of collocation yields equations for the Fourier
amplitudes, a gigantic, dense linear system. On the NERSC computer SEABORG, a 3,000processor IBM SP, a typical run requires 8 hours of processor time at approximately 1.7
teraflops, and memory of 750Mbytes/processor or 1,200 Gbytes total. Inputs are the MHD
equilibrium quantities (magnetic field geometry, plasma pressure in space), profiles of density
and temperature of a large number of particle species, distribution functions for non-Maxwellian
species (slowing-down distributions or superpositions of Maxwellians), and wave parameters,
descriptions of antenna geometry etc. The outputs are Fourier amplitudes of wave fields andplasma response obtained by post-processing: profiles of power deposition, driven current etc.,
obtained from the moments of f (evolution of f 0 + f 2 ). Calculating wave sources is in itself atask of integrated modeling, as this background evolution, on the same time scale as transport
processes, in turn affects the wave absorption. This full integration in calculating wave sources is
not presently achieved. A SciDAC project is making progress on integration, but the scope is far
short of what is needed
E. SUMMARY: THE ROLE OF SOURCE MODELS IN COMPREHENSIVEINTEGRATED SIMULATION

For whole device modeling, source models covering a range in levels of description will beessential elements of any comprehensive device model. We need a 3-D solution of the drift
kinetic equation including RF effects, radial transport, and finite orbit width effects.
In focused integration initiatives, source modeling will be essential. A few examples are:
• Particle source modeling is already a key element of edge studies. There are strong
interactions between the edge plasma and RF launching structures – RF fields perturb
the edge plasma, and edge plasma characteristics influence the wave spectrum launched.
This requires integration of RF models (particularly the antenna model) with an arsenalof edge modeling codes.
• In island formation, heating by RF waves at rational surfaces and energetic particles
produced during RF heating can stabilize island growth. RF driven currents can either
stabilize or destabilize island growth. This requires integration of RF models with MHD
and transport near the island (See Sec. VII).
FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
37•

To understand the role of radial electric fields on transport and flows in stellarators, weneed 4-D or 5-D solution of the drift kinetic equation including RF effects.

VIII. FUSION SIMULATION PRESENT CAPABILITIES STATUS
T his section summarizes the current status of integrated computational modeling andsimulation of toroidal confinement fusion systems. The intent in this section in the present
document is to provide a general perspective of the status of this very active and mature field.
This section responds to the explicit request in the subcommittee’s charge letter to report on thestatus of fusion simulation capabilities, and is the context from which the FSP is defined. The
present fusion simulation capabilities form a significant part of the critical underpinning of the
FSP.
There are over 50 major toroidal physics design and analysis codes being maintained by themagnetic fusion community. The major multi-user codes are depicted in Figure A11, whichshows how they divide into groups, and indicates with arrows the flow of information from one
code group to another.
The axisymmetric free boundary equilibrium codes solve the force balance equation by
calculating the poloidal magnetic flux in cylindrical coordinates for given pressure and current
profile parameterizations. These can be used to define the boundary for the inverse equilibrium
codes. There are also two major fully 3-D equilibrium codes in use: VMEC, which assumes theexistence of good magnetic surfaces a priori, and works in a coordinate system based on these,
and PIES, which calculates the existence of surfaces as part of the solution, if they exist.
The collection of linear macroscopic stability codes maintained by the MFE community is quite
mature and can assess the stability properties of a given equilibrium with respect to both ideal
and non-ideal (resistive) MHD, including the effects of an energetic particle component.

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
38

Major U.S. Toroidal Physics Design and Analysis Codes
Nonlinear Gyrokinetic
Continuum
globalGyrofluidVacuum &Conductors

Flux tube

VACUUM

Flux tube

GYRO

GS2

GTC

GRYFFIN

VALEN

Free Boundary
Equilibrium

Flux tube
SUMMITEFIT

Non-Ideal
high-nMHD- + particles

low-n

low-n

high-nPEST-I,II BALLOON PEST-III

NOVA-K

HINST

MARS

ORBIT

NOVA

CAMINO

DCON

BALOO

GATO

intermediate-n

Linear high-n
gyrokinetic

TEQ

FP-CodeAntenna

Ideal
low-n

Particle-in-cell
globalLinear Stability

TSCEQ

CQL3D

FULL

ELITE

RANT3DRF Heating
& CD
AORSA

TORCH

ORION

LSC

TORIC
METS

TORAY
CURRAY

2D transport
TSC

InverseEquilibrium 3D Nonlinear MHD

TRANSPJSOLVER

WHIST

TOQ

ONETWOESC

CORSICA2D plasmaneutrals

B2

DEGAS

NIMRODUEDGE

EIRENEFAR

3D plasma

Static Time -Dependent
PIES
M3D
VMEC

Plasma EdgeVMEC2DBOUT

POLAR2DBALDURdenotes parallel MPI code

Fig A11. Major U.S. toroidal physics design and analysis codes used by the plasma physics community.

The nonlinear codes fall into four major groupings. In descending order of the frequencies
addressed, these are the: 1) RF Heating and Current Drive codes, 2) the Nonlinear Gyrokinetic
codes, 3) 3-D Nonlinear Extended MHD codes, and 4) the 2-D Transport codes.
The RF Heating and Current Drive codes calculate the propagation of electromagnetic waves ofa given frequency through prescribed background plasma, including reflection and absorption.
The codes are of two major types: ray-tracing (or geometrical optics), and full wave (global
solution). There are also depicted antenna and Fokker-Planck codes, which are closely coupledwith the RF codes, and provide boundary conditions and background distribution functions.
The RF codes are designed to calculate the instantaneous heating and current-drive profiles for agiven plasma equilibrium subject to a given RF oscillator source and antenna.
The gyrokinetic codes are based on an analytic reduction of the full 6-dimensional plus time
plasma distribution function obtained by averaging over the rapid gyro-motion of ions in astrong magnetic field, and by neglecting the displacement current in the Maxwell equations toremove “light waves” from the system. These codes are appropriate for studying 3-D turbulent
transport in a background system with fixed profiles.

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
39

The 3-D nonlinear Extended MHD codes are based on taking velocity moments of the
Boltzmann equation to yield 3-D magneto-fluid equations for the evolution of the averageplasma velocity, density, and pressure, along with a closure procedure. These codes areappropriate for describing global stability phenomena such as sawteeth oscillations, magnetic
island evolution, and plasma disruptions.
The 2-D transport codes presently form the core capability in our community for integrated
modeling. There are six major codes, with considerable overlap, that exist largely for historicalreasons. These codes are all based on the Grad-Hogan evolving equilibrium description where
the inertial terms in the momentum equation are neglected and the remaining MHD equations
are averaged over the flux surfaces, where they exist.
The 2-D transport codes are all very modular. They are each a collection of equilibrium
modules, transport modules, solvers, and source and sink modules representing Neutral Beam
Injection (NBI) and RF heating, pellet and gas injection, impurity radiation, and the effects ofsaturated MHD activity such as sawteeth and islands. These codes have recently benefited fromthe National Transport Code Collaboration (NTCC), which has formed a modules library sothat modules taken from individual codes can be exchanged and shared (see, e.g.,
w3.pppl.gov/NTCC for more details on this). While these codes address integrated modeling,
the individual modules represent simplified reduced descriptions of the full three-dimensional
physical phenomena being modeled.
A summary of these four major code groups and the timescales being addressed by them has
been given in Fig. A6. The RF codes address frequencies of order the ion cyclotron frequency,Ωci , and above, up to the electron cyclotron frequency Ωce. (Again, recall the interchangeable
notation for these two quantities, ω ci and ω ce .) The gyrokinetics codes typically take time steps
about 10 times longer than the ion cyclotron frequency, whose motion is analytically averaged
over, and are normally run for 103 to 104 time steps to calculate stationary turbulent fluctuationlevels. Recent additions to these codes to include some electron timescale phenomena bring in
the electron transit time, which lowers the maximum time step. The Extended MHD codes needto resolve phenomena occurring on the Alfvén transit time, τA, although most codes are at least
partially implicit to avoid a strict restriction on the time step based on this. These codes cannormally run 104 to 105 time steps to address MHD phenomena such as sawteeth and island
growth. The 2-D transport codes are very efficient and can take many long time steps to model
the entire discharge. However, the 2-D edge transport codes need to resolve the parallel
dynamics and use a time-step based on the parallel sound-wave propagation near the edge region.
The calculations now being performed with the gyrokinetics codes, the Extended MHD codes,
and the RF codes, are straining the limits of the existing computing capabilities and capacities.
For example, recent attempts by the core-turbulence gyrokinetics codes to include both electron
and ion dynamics in a self-consistent simulation require upwards of 104 processor-hours (overone processor-year) on the IBM SP3 at NERSC to generate one result for a set of fixedbackground profiles. Similar times are required by the Extended MHD codes to calculate thegrowth and self-consistent saturation of a neoclassical tearing mode. Thus, we can take solace inthe fact that the capability is mostly in place, but must deal with the fact that the computational
requirements for a fully integrated 3-D comprehensive simulation capability are truly daunting.

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
40

Examples of fundamentally important experimental phenomena that involve 3-D physicalprocesses that cross theoretical boundaries and thus cannot adequately be addressed by thepresent suite of above-described codes include:

Pedestal physics – A description of the transport barrier that forms in the region of the plasma
between the core and the edge, and of the associated edge localized relaxation events;
Long time scale profile evolution – A way to self-consistently evolve the global profiles of plasma
temperature and density on the energy-confinement time scale from turbulent transport and in
the presence of magnetic islands and other MHD phenomena.
Edge transport: A description of long-mean-free-path particle and heat transport outside theclosed magnetic flux surfaces, on the open field lines that impact the first wall or divertor andinvolve multi-phase physics
Self-consistent heating and current drive: A fundamental model of the interaction of Radio
Frequency (RF) waves with plasma in the presence of plasma turbulence.
Sawtooth phenomena – Internal MHD-type modes in the hot core of tokamak plasmas forwhich fast ion and kinetic effects are clearly relevant experimentally but are only beginning to beaddressed computationally.
Island physics: incorporation of the effect of 3-D island formation on equilibrium evolution andturbulent transport
These are but some of the important problems to be addressed by the integrated simulation
initiative.

IX. FOCUSED INTEGRATION INITIATIVES – CHALLENGES ANDOPPORTUNITIEST he following Subsections A-D expand upon the descriptions of the candidate focused
Integration Initiatives (FIIs) described in the main text. These are based on reports from working
groups at the September 23-24, 2002 Community Workshop, each concentrating on a particular
FII.
A. PLASMA EDGE
The boundary or edge-plasma of fusion devices plays a vital role in their operation. The edgeplasma region is generally considered to be the region where substantial two-dimensional (2-D)
or 3-D variations can occur in the plasma, neutral particle, and magnetic equilibrium quantities.
In addition, owing to the lower plasma temperature and proximity to material surfaces: neutralgas, sputtered impurities, and atomic line-radiation can become important components. Thus, arich variety of potential interactions can take place in this region.
The ISOFS edge discussion group identified four elements that are key to successful operation ofan MFE fusion device:(1) predicting conditions and properties of the pedestal energy transport barrier just inside
the magnetic separatrix;
FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
41

(2) understanding plasma/wall interactions for particle recycling and wall lifetime from highenergy fluxes;(3) controlling tritium inventory including co-deposition; and(4) controlling wall impurity production and transport into the plasma.
A number of models of varying sophistication exist to describe these processes. Some models
already provide a level of coupling, e.g., hydrogen transport, recycling neutrals, impurity
sputtering, and impurity transport codes. However, many of the constituent models need
improvement, and much more inclusive couplings are required to self-consistently predict theedge-plasma behavior.
In setting priorities for the edge region, it was felt that the short-term (five-year timescale) work
should focus on a better understanding of what controls the suppression of plasma turbulence toproduce a transport barrier in the pedestal region (#1 above), and its associated impact on plasma
profiles. The ability to predict the behavior of the edge pedestal barrier is essential for projecting
the net fusion output of MFE devices as discussed extensively at the Snowmass 2002 Fusion
Study Workshop. Presently, the key parameter that is believed to control the core fusion output
is the plasma temperature at the top of the pedestal; this parameter is now either extrapolated
from existing experiments or assumed.
The development of a predictive model of the pedestal requires two key advances. The first is theinclusion of kinetic effects in the 3-D turbulence simulations such that the physics model can
span the region from the nearly collisionless portion at the core side of the pedestal to the more
collisional region near and outside the magnetic separatrix. Presently, turbulence codes fall into
two classes: one class is the collisional fluid codes allowing three spatial dimensions and including
the strong magnetic shear region in the edge; the other class is the kinetic particle or continuum
codes using three spatial and two velocity-space dimensions to capture the physics of plasma withlong mean-free paths, but allowing only very simple plasma geometry. These codes are very
compute-intensive, requiring parallel computers with run times of many hours or days (> 2000
node-hours). The computational challenge here is to obtain a generalized code or coupling ofcodes that can treat the full range of collisionality through the pedestal for geometries relevant to
the plasma edge.
Another component needed for the pedestal is the description of periodic edge profile relaxations
called Edge Localized Modes (ELMs), which are believed to be caused by long wavelength MHD
instabilities. The transport barrier is intermittently interrupted by these modes. The peak heat
flux to material surfaces from such events is a serious concern for a reactor-sized device. Inclusionof such modes in the microturbulence codes should be considered, but requires a larger spatial
mesh and an extension of the physical model to capture this MHD phenomenon. An alternative
approach is to use a separate code for the turbulence from MHD modes, or to use a linear
stability code to enforce marginal stability of the plasma pressure profile for long-time transport
evolution.
A second key advance thought possible in the five-year time frame is to couple the turbulent
plasma fluxes with a transport calculation of the profile evolution. Such coupling is the heart of apredictive model. This is the plasma-edge analogue of the “Turbulence on the Transport
Timescale” FII described in the next section. The turbulence drives plasma fluxes that, together

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
42

with sources such as neutral ionization, largely determine the plasma profiles; since these plasma
gradients provide the free energy to drive the turbulence, it is essential that the turbulence andprofile evolution be coupled together. Plasma-wall interactions of sputtering and recycling canmake the profile evolution time substantially longer than the turbulence saturation time.
Devising stable and efficient coupling schemes is a major task here, which may have spin-off torelated couplings for other regions of the fusion plasma and beyond. Experiments and turbulencesimulations often show large, intermittent transport events that need to be properly
characterized.
Beyond the five-year time frame, edge-plasma modeling should work to develop and couple
more detailed models of the plasma-wall interactions. These include models for sputtering,
recycling and re-deposition, and transport of neutral gas impurities and radiation. In addition,
the edge model should look inward to couple with core physics inside the pedestal region. A
vision is that a successful edge model could eventually incorporate the core region models for awhole-device simulation.
B. TURBULENCE ON TRANSPORT TIMESCALES
The nature of the problem to be considered can be summarized as follows: The “anomalous”
transport of mass, energy, and angular momentum in toroidal MFE devices is dominated byfluxes driven by plasma turbulence. There is a significant disparity of scales, especially timescales.
This is a highly coupled system: the plasma pressure and densities are very nearly constant onmagnetic flux surfaces, and can thus be considered to be one dimensional (1-D) functions of theflux coordinate ψ; one-dimensional transport equations can be constructed to describe theevolution of these profiles, whose gradients drive turbulence, while the turbulence (3-D,
anisotropic) produces the fluxes that drive the evolution of these 1-D profiles.
The object of this Focused Integration Initiative (FII), is to bridge the range of temporal andspatial scales so as to compute the above coupled system self-consistently, as opposed to justcomputing 3-D fine-scale turbulence with fixed background profiles, or computing 1-D
transport with highly reduced theoretical or empirical models of the turbulent fluxes.
The single overarching science issue and goal is the self-consistent calculation of global
confinement from first-principles physics. The initial (easier) focus would be to determine
steady-state confinement; following time evolution on the transport timescale is conceptually nomore difficult but is more computationally demanding. The achievement of the steady-state goal
would, as a side benefit, enable optimization studies. Important issues like simulation of both
steady-state and time-dependent versions of internal transport barriers are subsets of the overall
goal.
There are a number of candidate approaches for achieving the goal. These include:(1) direct coupling of the transport and turbulence equations;(2) an expert system working together with a smart database of past turbulence simulation
results and/or analytic models for transport coefficients, which detects when transport
has moved outside of the domain of applicability of the database and triggers new
turbulence simulations to extend the domain of applicability;

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
43

(3) developing a viable gyrofluid closure and directly integrating the resulting fluid
equations; and(4) projecting long-timescale evolution from moments of the extrapolation of particle
weights in a gyrokinetic particle-in-cell code.
This candidate FII entails many disciplines in physics and computation. It also entails a numberof generic integration issues, including: coupling of disparate timescales; coupling of descriptions
with different dimensionality; interfacing to experiments for validation; the need for a flexible
framework/environment to support experimentation with multiple approaches; visualization;
expert systems to detect failure of an integration component; software standards; and databasemanagement. All of these issues present opportunities for collaboration between thetheoretical/computational plasma physics and the computer science/math communities.
Such an initiative would also need a strong theory component. Most importantly, there needs tobe a parallel effort in analytic turbulence theory to provide crosschecks and understanding.
Theory will also be needed for the gyrofluid approach and other formalism extensions, and fordevelopment of synthetic diagnostics.
There are clearly opportunities for further linkages and progress toward a full integrated
simulation. Turbulence and MHD dynamics may need to talk to one another directly as
opposed to via the intermediary of a transport code, but this coupling may be a logical extensionof the transport-turbulence coupling activity, as the necessary flows of information are similar.
Another logical extension is coupling of the turbulence equations to real sources. But, perhaps
most importantly, if the transport-timescale part of transport-turbulence coupling is handled by
an existing or new (developed under the initiative) full-device integrated modeling/transport
code, this activity immediately has access to many other developments in whole-device
integration.
C. GLOBAL STABILITY
Global stability issues play a central role in determining the optimal operating regime of fusion
devices, and in describing their time evolution. It is well known that under some operating
conditions, an experimental discharge can spontaneously transform from a symmetrical stablesystem exhibiting good confinement into one that exhibits symmetry-breaking oscillations andpoor confinement or becomes unstable and disruptive. These events are known as sawtoothoscillations, Neoclassical Tearing Modes (NTM), disruptions, Edge Localized Modes (ELMs), or
Resistive Wall Modes (RWM), among others. The predictive calculation of the onset andevolution of these events is the over-arching scientific question in this FII.
The long time evolution of tokamak discharges is often described by nested magnetic flux
surfaces and two-dimensional axisymmetric force balance coupled with diffusive transport of the
one dimensional surface averaged thermodynamic and field quantities. The transport coefficientsused in the diffusive model are typically analytic approximations to the results of threedimensional, sub-grid-scale turbulence that is described by kinetic theory. The properties of the
turbulence are determined by the axisymmetric plasma profiles, which in turn are affected by thetransport coefficients. The first-principles coupling of the transport and kinetic turbulence

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
44

models is a formidable problem requiring integrated modeling as described in the previous
section.
However, the force balance assumed in the transport model can occasionally become unstable
and yield dynamics that evolve on an intermediate time scale that is much shorter than thetransport time scale, but much longer than the time scale of the kinetic turbulence. Thesemotions result in three-dimensional magnetic perturbations that break both the underlying
axisymmetry of the overall configuration and the nested topology of the magnetic field. Their
nonlinear evolution can strongly affect the confinement properties of the magnetoplasma system,
and can re-arrange the plasma profiles on time scales faster than that described by diffusivetransport.
At relatively low temperatures, these dynamics are well described by r e s i s t i v emagnetohydrodynamics (MHD), a mathematical model in which the plasma is treated as anelectrically conducting fluid subject to electromagnetic body forces. Solutions of this model arecomplicated by a wide separation of space and time scales, and by the inherent high degree ofanisotropy that occurs in a magnetized plasma. At the higher temperatures that occur in modern
tokamaks, kinetic effects both parallel and perpendicular to the magnetic field introduce
important physical processes that can affect the global MHD evolution of the plasma. The
challenge is to develop mathematical and computational models that include these kinetic effects
while retaining the computational tractability of the fluid model. Such models are collectively
called extended MHD.
There have been two approaches to developing extended MHD models, as described in Section
IV. One is to take analytic moments of the higher dimensional drift kinetic equation to obtain
expressions for the higher order fluid closures (e.g., the heat flux and the components of the stresstensor) that capture the kinetic effects within the fluid model. The other is to solve the kinetic
equation by sub-cycling within a time step of the MHD model. Velocity moments of theresulting distribution function can be taken numerically to calculate the required closures. While
this approach may be sufficient for minority (or possibly even majority) ion species, it difficult toenvision it being applied to electron dynamics because of their small mass and rapid dynamics.
There it is likely that we must continue to rely on analytic closures.
A full modeling of kinetic effects on MHD evolution, and capturing the effects of the MHD
relaxation of the profiles within longer time scale transport calculations, requires an integrated
simulation. For example, consider the long time scale modeling of a burning plasma
configuration. The axisymmetric profiles evolve in response to the turbulent transport and thesources of mass, momentum, and energy from RF antennas and neutral beams, which are in turnaffected by the evolution of the profiles. Throughout the simulation the stability of theconfiguration can be monitored. The stability is be affected by the profiles and the presence ofthe energetic alpha particles that are produced by the fusion reactions. When the configuration
becomes unstable it can be used as the initial conditions for an extended MHD simulation. Thismodel would predict the dynamics of the three-dimensional magnetic perturbations anddetermine their effect on the underlying profiles. In a burning plasma, this calculation would
require integrated models for majority ion and electron damping, and for minority (alphaparticle) dynamics. The resulting modified profiles could then be returned to the transport
calculation for further evolution.

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
45

To make progress on this Focused Integration Initiative would require the talents of theoreticaland computational plasmas physicists, applied mathematicians, and computer scientists. The
existing Extended MHD codes need to be enhanced to include higher resolution and a more
rigorous mathematical framework for coupling MHD events with both microscopic effects andwith long timescale evolution of profiles, eventually including such effects as plasma rotation andedge effects. The applied math community will be called upon to provide novel methods to deal
with the time and spatial scale separation and with the extreme anisotropy; for example moving
or adaptive grids, high-order or spectral elements, and nonlinear equation solvers. Since there isno agreed upon “best” method for solution of the Extended MHD equations or for the
couplings, there is a need for a computational framework that allows rapid prototyping. Many ofthe problems faced by this initiative are generic to the larger integration problem. Indeed, since
MHD phenomena are intermediate between very fast and very slow phenomena, this FII is a
microcosm of the entire project.
D. WHOLE DEVICE MODELINGThe distinguishing feature of the Whole-Device Modeling (WDM) FII is that from the outset it
will provide a model of the entire device for the whole discharge timescale. Because of this scope,
many of the models of individual systems are necessarily very simple. It is envisioned that these
simple models will be capable of being replaced by more complete and accurate models as they
become available and as is warranted by the application.
The state-of-the-art of whole-device complete-shot modeling at present is represented by thenational array of 1-1/2-D transport codes described in Secs. VI and VIII. Those codes have
many features that would be required for a final product WDM. They employ a formal
separation of time-scales between the rapid (Alfvén) time on which 2-D magnetic equilibria areestablished, and the much slower time on which heat, particles, and angular momentum, aretransported as 1-D surface functions across the magnetic surfaces. They also incorporate many
features of a WDM: a hierarchy of models to describe particular aspects of physics, with tradeoffs between speed and accuracy; connection to experimental databases; predictive andinterpretive modes. Some are also equipped with sophisticated user interfaces.
There are at least three distinct thrusts in a WDM FII. The first is to extend the accuracy andreach of the physics modules available in the existing codes; the second is to extend these codes
to fully 3-D geometry so they are applicable to stellarators; and the third is the development of asuitable modern computing framework architecture that would allow this effort to couple to thefruits of the other FIIs for the final WDM code.
The first of these involves the development and application of algorithms for coupling “bestphysics” modules to the surface-averaged long-timescale transport equations described in Sec. VI.
The scientific goals of such couplings are many and diverse: the accurate prediction of heating,fueling, current drive, and confinement on the transport timescale; the ability to accountaccurately for sawteeth, tearing, and wall modes, in order to assess, avoid, or control them;increased understanding of the complex interaction between the plasma edge and coreconfinement; and the ability to model and develop machine feedback and control systems withina simulation code environment. Many of these capabilities exist already at some level, and thisFESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
46

activity can be viewed as extending these and taking them to the next level. Extensions would
include incorporating rotation and non-Maxwellian distributions into the equilibrium, and selfconsistently incorporating multiple heating and current-drive systems into the formulation.
The second major thrust involves coupling 3-D plasma equilibria and device geometry to 1-D
surface-averaged transport and the many source and boundary modules that make up a Whole
Device Modeling code. This can be viewed as an extension of the present 2-D equilibria + 1-D
surface-averaged transport codes, but it is a non-trivial extension. The calculation of 3-D
equilibrium is itself a research topic, and issues of existence of surfaces, magnetic islands, andfield line stochasticity must be dealt with on both a fundamental mathematical and a practicalcomputational level. It will also be necessary to extend current methods to include plasma
rotation and other non-ideal effects. Developing a system to define the complex 3-D input
geometry and to manage the extensive output is itself a challenge.
The third major goal of a WDM FII is the design of the code architecture for the final WDM
product. This is obviously critical in many respects. Although it may not speed science results
within the five-year timeframe of the FIIs, it is important for the success of the global initiativethat the framework itself gets a timely start. This thrust, in particular, will rely heavily on
computer scientists to develop an extensible framework that meets the many needs of the
projects.
By its nature, WDM is an integrated activity, and initiatives in this area will overlap other FIIs.
While development and testing of particular couplings will be carried out as narrowly focused
projects in the other FIIs, integration into a WDM code is a practical requirement for a close
connection to experiment, or for an ability to simulate a proposed new machine on transport
timescales. It is also possible for a WDM code to serve as the 1-D transport solver throughout
the development of any of the new couplings. This suggests coordinating the WDM FII closelywith the turbulence, the MHD, and the edge FIIs.

X. VALIDATION REQUIREMENTS IN FUSION TOPICAL AREAS
W hat follows are examples of phenomena that one could anticipate being important forvalidation tests of fusion simulation codes in various topical areas, with particular regard topossible Focused Integration new capabilities. This is not a statement of the outstanding physics
issues - but a list of some calculable and measurable (at least in principle) quantities that areimportant for testing. They are described beginning with the most general and global to the mostdetailed and local and most stringent. This ordering also reflects the difficulty in making suchmeasurements.
A. TRANSPORTThe coarsest level of agreement between a transport simulation and an experiment would be tomatch the total stored energy given the input power and other machine parameters. This is only
an interesting comparison if the model has few if any free parameters that have been calibrated
against existing databases. One should not expect the scaling laws themselves to be derivable as
FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
47

the engineering parameters used are only proxies for the relevant physics variables. A higher levelof agreement could be assessed by comparison of profiles given fluxes or models for the sources.
At a minimum, one is interested in the width and height (or gradient) of the edge pedestal andcore temperature gradients. Eventually this comparison must be extended to all transportchannels, ion and electron thermal, particle (including impurities) and momentum transport.
Recent work by ITER expert groups have shown however that a wide range of codes, based on
physical, semi-empirical, or wholly empirical models can achieve essentially identical
performance in matching temperature profiles for a moderately large set of discharges from avariety of machines. This suggests that agreement at this level is insufficient to validate aparticular code. More challenging are comparisons of transient behavior including thresholds
and dynamics of transport barriers. Many widely observed features of transient transport cannot
be explained by the current generation of models. Ultimately, transport models must be
validated at the level of turbulence and turbulence dynamics - comparisons that will be paced by
the development and deployment of fluctuation diagnostics and analysis techniques. Adding to
the difficulties is the prediction that fluctuations vary significantly over the poloidal cross section.
Synthetic diagnostics, numerical analogies to the experimental diagnostics need to be developed
and adapted to simulation code outputs.
B. MHDThe first test of MHD codes is their ability to reproduce experimental stability limits oroperational boundaries for ideal and resistive instabilities, including kinks and ballooning modes,
edge localized modes (ELMs) and other edge relaxation phenomena, and resistive wall modes.
Especially interesting are modification to the stability boundaries from non-ideal or non-linear
effects. A second parameter for comparison is the growth rates of unstable modes, first in theirlinear phase then in the non-linear phases including the calculation of the saturation levels forthese modes. More generally one can look for agreement with the computed eigenvectors andeigenvalues including non-linear mode structures. For example, one should be able to predictvarious ELM types and their non-linear extent. Large-scale dynamics of disruptions can becompared including halo and eddy current distributions in real machine geometry, runaway
populations and so forth. Dedicated experiments should also test explicitly the extensions toideal MHD including neoclassical, two-fluid, flow, finite Larmor radius and other kinetic effects
and non-linear interactions with profiles and transport.
C. RADIO-FREQUENCY HEATING AND CURRENT DRIVE
Validation of RF models is challenging, as the important quantities are particularly difficult tomeasure. Testing begins with global quantities like overall heating and current drive efficiencies
and proceeds to comparisons with deposition profiles for heat, current, and flow velocity.
Beyond these measures, an essential element is the verification of RF waves inside the plasmas.
To the extent possible, the two-dimensional fields of wave amplitude and wave number shouldbe measured and compared to code predictions. The position of mode conversion layers shouldbe verified along with the propagation of the outgoing waves. Antenna/edge-plasma interactions
and the influence of plasma fluctuations on launched waves will challenge both simulations andmeasurements. Finally, one will need to compare codes with experiments that test the models ofwave-particle interactions. For this, one would need to measure velocity-space distributions alongwith the wave fields and plasma profiles.

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
48

D. EDGE/SCRAPE-OFF-LAYER/DIVERTOR
Edge modeling involves a wide range of physical effects including those common with the corelike transport and MHD as well as neutral dynamics, atomic processes and plasma-materialinteractions. Fortunately the diagnostic challenge is not quite as severe in the edge plasma as it isin the core. Assessment of edge transport models is similar to that in the core with comparisonsinvolving profile and fluctuation measurements, though with perhaps greater emphasis on
poloidal variation. The energy source from the core to the edge is expected to vary strongly
around the poloidal circumference leading to predictions of flows and potentials that must be
verified. As in the core, profiles from energy, particle, impurity, and momentum transportshould all be compared along with the appearance of self-generated flows. Unlike the core, thesources and sinks for particles, including impurities, are not well characterized. The position andprocesses for these sources needs to be included in models and tested by experiments. Depositionof impurities and co-deposition of hydrogenic species is a particularly important issue forbenchmarking. The role of neutrals is complex, affecting all transport channels through bothclassical and turbulent processes. The computed three-dimensional distribution of neutrals
should be verified and the role of various neutral transport mechanisms tested. A basic model for
fueling including sources, neutral dynamics, particle transport needs to be developed andcompared with experiments. Finally, important interactions with MHD physics in phenomena
like ELMs and the tokamak density limit provide particularly stringent tests for integrated
models.

XI. BACKGROUND – ORBITS AND INSTABILITY MECHANISMS
H ere, we explain some of the basics of plasma motion embodied in Eq. (A-1). For a more
complete explanation, one should consult any plasma physics textbook. The basic equation,
describing the motion of a charged particle in electric and magnetic fields, reflecting the force
term in Eq. (A-1) ismdv= q(E + v × B) ,dt

(A-28)

where q is the particle charge. The basic motion in a uniform magnetic field is gyration about thefield line with frequency ω c = qB m and radius ρ c = v ω c . Additional forces and non-uniformities
produce drifts in a direction perpendicular to both the force and the magnetic field. Thus, an electricfield E or a gradient ∇B perpendicular to B reduces the orbital radius on one-half cycle and increases
it on the other, resulting inv ⊥2 B × ∇B
E×B
, and v d =,vE =2ω c B 2
B2(A-29)

the E×Β drift and the grad B drift, respectively. If the magnetic field is curved, a similar
expression is obtained for the curvature drift, which is proportional to v||2 . In addition to theseparticle drifts there are fluid drifts in confined plasmas (i.e., plasmas with radial density and
FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
49

temperature gradients), owing to the fact that at any one radius, there are more particles moving
in one direction perpendicular to the gradient than in the reverse direction. This yields a net
fluid velocity known as the diamagnetic driftv*a =ρ a v a B × ∇n a,
Bn a

(A-30)

where v a is the thermal velocity of species a. A similar drift proportional to the temperature
gradient is also important in the kinetic theory. For the orderings employed in plasmacomputation, an important feature of all these drifts is that for thermal particles they are smaller
than the thermal velocity by a factor ρ a L << 1, where L is a macroscopic length, e.g. of orderρ * . All plasma waves, instabilities, and transport phenomena are related in some way to thesedrifts, which are illustrated in Fig. A12.

Fig. A12. Ion drifts in electric and magnetic fields.

The existence of waves and instabilities is closely related to the fact that the grad B drift and the
diamagnetic drift have opposite signs for ions and electrons (we are taking ρ a and ω ca to havethe sign of qa ), while the E×B drift is the same for each species. For example, Fig. A13 provides
a simple picture of instability in bad curvature. If a perturbation forms in a region of highpressure gradient, then curvature and grad B drifts cause a charge separation. The resulting
electric field drives E×B drifts that amplify the perturbation.

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
50

Fig. A13. Bad curvature instability

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
51

XI. GLOSSARY
Alfvén wave: A plasma wave, which involves bending or compression of the magnetic field. The
Alfvén time scale is the interval for an Alfvén wave to traverse the plasma. Magnetosonic waves
are a type of Alfvén wave.
Alpha heating: In a fusion power plant, energetic alpha particles and neutrons are created by thefusing of deuterium and tritium nuclei. As a charged particle, the alpha particle is unable to crossthe confining magnetic field and gives up its energy as heat to the plasma.
Ambipolarity: Mass transport in plasmas is ambipolar, that is the fluxes of electrons and ions are
virtually the same. The densities of negatively and positively charged particles that compose aplasma are almost in perfect balance, leaving the plasma essentially neutral (termed quasineutrality).
Aspect ratio: In a toroidal device, the ratio of the major radius to the minor radius.
Avalanche: A sudden macroscopic event in which energy or particles can be distributed acrossthe medium. Often described by simple mathematical models.
Bad curvature: See Curvature.
Ballooning: A local instability, which can develop in the tokamak when the plasma pressure
exceeds a critical value. It is analogous to the unstable bulge that develops on an over-inflated
pneumatic inner tube.
Banana: The shape of a trapped particle orbit (banana orbit) projected onto a poloidal plane.
The shape results from drifts away from a magnetic surface.
Beta: The ratio of plasma pressure to magnetic pressure. An essential dimensionless parameterfor magnetized plasmas. Denoted by the symbol β.
Boltzmann equation: An equation of motion in phase space (position and velocity). The
Boltzmann equation (and its collisionless version, the Vlasov equation), is the starting point formuch of plasma physics and is the fundamental equation for kinetic theory.
Bootstrap current: Currents driven by collisional transport effects in toroidal plasmas (seeneoclassical).
Braginskii equations: Plasma fluid equations separately describing electron and ion motion.
Fluid equations are obtained by calculating velocity moments of the Boltzmann equation andspecifying a closure condition.
Closure: Mathematical scheme by which a hierarchy of equations is truncated. This usually
involves expressing higher moments in velocity space in terms of lower moments.
C-Mod: A compact high field tokamak experiment.
FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
52

Confinement: Property of magnetic fields in preventing loss of energy and particles; degree ormeasure of this property, as in “confinement time,” for example.
Curvature: In toroidal configurations, field lines are necessarily curved, giving rise to particle
drifts. In bad curvature or unfavorable curvature regions, these drifts give rise to instability (seeFig. A13). In good curvature regions, the drifts are stabilizing, owing to the reversal of the relative
direction of curvature and pressure gradient. The grad B drift acts in a similar way.
Cyclotron frequency: Same as gyrofrequency.
Detached plasma: Cold dense plasma in the divertor chamber, separated from the walls by
neutral gas. This is characteristic of a particular mode of divertor operation.
Diamagnetic drifts, waves, etc.: Refers to plasma dynamics driven by the effect of a plasmapressure gradient in a magnetic field.
DIII-D: A medium scale, strongly shaped tokamak experiment.
Disruption: The abrupt termination of a tokamak plasma through the growth of large
amplitude MHD instabilities. Control of disruptions is a critical problem for this type ofconfinement device.
Divertor: Region outside the plasma core with open field lines leading to a chamber somedistance from the plasma.
Drift: The motion of a particle or fluid when subjected to a force or gradient perpendicular to
the magnetic field. The drift is perpendicular to both the field and the other force. Examples are
curvature drift, grad-B drift, E×B drift, and diamagnetic drift. See Appendix A8.
Drift wave: A type of plasma wave arising from the presence of density and temperature
gradients across magnetic field lines. Turbulence of various types of drift waves is believed to be
responsible for anomalous transport in toroidal plasma experiments. ITG and ETG modes areexamples of drift waves.
ECRF: Electron Cyclotron Range of Frequencies – Refers to radio frequency heating and current
drive using waves close to the electron cyclotron frequency.
ELM: Edge Localized Mode – a phenomenon that relaxes the pressure gradient of a plasma
device operating in H or High confinement mode.
Equilibrium: Usually refers to MHD equilibrium – a steady state solution of the MHD
equations. While confined plasmas may be approximately in local thermal equilibrium, they are
not in global thermodynamic equilibrium.
ETG: Electron Thermal Gradient (modes) – A type of “micro” plasma instability that may be
responsible for turbulent transport by very short (electron cyclotron radius) fluctuations.

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
53

Field lines and flux surfaces: Imaginary lines marking the direction of a force field. These map
out surfaces (to which plasma particles are approximately constrained) called flux surfaces.
Finite Larmor radius: FLR. A mathematical expansion or approximation in which terms oforder the Larmor radius (gyroradius) divided by a macroscopic scale are retained.
FIRE: Fusion Ignition Research Experiment – A proposed burning plasma experiment, which
would operate with very high magnetic fields.
FLR: See Finite Larmor radius.
Flux: See Magnetic flux, Transport flux.
Good curvature: See Curvature.
Grad-Shafranov equation: A steady-state MHD equation whose solution yields axisymmetricequilibria.
Gyrofluid: A set of fluid equations that treat ions and electrons separately and retain certain
kinetic effects through their closure conditions.
Gyrofrequency: the frequency at which a particle gyrates around a magnetic field line; also
called cyclotron frequency.
Gyrokinetic: Kinetic equations derived from the Boltzmann equation, which is averaged overthe fast cyclotron (gyro) motion of the plasma particles reducing the number of spatial
dimensions from 6 to 5.
Gyroradius: the radius at which a particle gyrates about the magnetic field; also called Larmor
radius.
Hall term: A term in the extended Ohm’s law in MHD, proportional to the crossed current andmagnetic field. The presence of this term, often neglected in simpler descriptions, leads towhistler waves and can strongly influence the rate of magnetic reconnection.
H-Mode: High (confinement) Mode – an experimental regime with a transport barrier orpedestal at the edge of the plasma.
IBW: See Ion Bernstein Wave.
ICRF: Ion Cyclotron Range of Frequencies - Refers to radio frequency heating, current and flow
drive using waves close to the ion cyclotron frequency.
Ion Bernstein Wave (IBW): – short wavelength electrostatic waves in the ion cyclotron range of
frequencies. Generated in the plasma by mode conversion.

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
54

Ion cyclotron wave: A plasma wave that propagates at frequencies above the ion cyclotron
frequency. See ICRF.
IPPA: Integrated Program Planning Activity – the current roadmap for the fusion energy
program.
Island – Magnetic Island: A three dimensional magnetic structure arising when closed flux
surfaces are perturbed by magnetic fields from plasma fluctuations or from errors in the vacuum
fields.
ITB: Internal Transport Barrier – a regime of strongly reduced transport over some part of thecore plasma.
ITER: International Thermonuclear Experimental Reactor – A burning plasma experiment
proposed as the next step in the international fusion program.
ITG: Ion Thermal Gradient (modes) – A type of “micro” plasma instability that may be
responsible for turbulent transport by medium scale (ion cyclotron radius) fluctuations.
JET: Joint European Torus – a large Tokamak experiment run jointly by a European
Community consortium.
Kinetic: Refers to aspects of plasma physics related to velocity distributions of electrons andions.
Kink Instability: A macroscopic (large scale) instability driven by plasma currents.
Larmor radius: Another term for the cyclotron orbit radius of electrons or ions in a magnetic
field.
LH: Lower Hybrid – RF waves of frequency intermediate between the electron and ion
gyrofrequencies and used for heating and current drive.
L-Mode: Low (confinement) Mode – the baseline for confinement in magnetic fusion devices,dominated by strong turbulent transport.
Lorentz force: The basic electromagnetic force on a particle, q(E + v × B) , where q and v are
the particle charge and velocity, and E and B are the electric and magnetic fields.
Lundquist number: The ratio of the resistive time scale to the Alfvén time scale. Plays a role for
MHD waves analogous to the Reynolds number in fluid dynamics.
Magnetic flux: The integral over a surface area of the magnetic field perpendicular to thesurface. Examples are the poloidal flux Ψ and the toroidal flux Φ.

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
55

Magnetic surface: A two-dimensional closed toroidal surface on which magnetic field lines lie.
In an axisymmetric device, such as a tokamak, the equilibrium consists of nested magnetic
surfaces.
Magnetosonic wave: Plasma sound wave, in which the magnetic field is compressed along withthe gas.
Major radius: Distance from the center line of the torus to the magnetic axis or center of theplasma cross section.
Maxwellian: Velocity distribution characteristic of local thermal equilibrium: f ~ exp(-mv2/2T)
where T is the local temperature. (See Equilibrium).
Mercier criterion: A mathematical criterion describing the local stability of MHD modes.
MHD: MagnetoHydroDynamics – A fluid description of magnetized plasmas.
Minor radius: distance from the center of the plasma cross section to the plasma edge.
Mode: 1) a plasma wave that has, at least approximately, constant wavenumber and frequency,
or can otherwise be characterized as a normal mode of some wave equation. 2) a mode ofoperation of a device, such as L-mode or H-mode.
Mode conversion: The process by which waves of one type are converted to another, occurringin a region where the frequency and wavelength of both waves match.
NBI: Neutral Beam Injection – A common method of plasma heating which employs intense
beams of neutral atoms.
NCSX: National Compact Stellarator eXperiment – a low aspect-ratio stellarator now in thedesign phase.
Neoclassical: The theory of collisional transport in toroidal geometry.
NTM: Neoclassical Tearing Modes – Resistive MHD modes driven unstable by bootstrap
currents.
Ohmic: Characteristic of dissipation due to plasma resistivity; Ohmic heating, in which acurrent is driven inductively in the plasma.
Pedestal: Region of steep gradients at plasma edge.
PIC: Particle In Cell – a numerical method for solving kinetic equations by following
representative particles moving in self-consistent fields. (Cell refers to the grid used in thecomputation.)
Poloidal: The short way around in a torus–denoted by θ (see Fig. III.1 of the main report).

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
56

Potato: The shape of a particle orbit (potato orbit) close to the magnetic axis – a distorted
banana, which see.
QPS: Quasi-Poloidal Stellarator: a low aspect-ratio stellarator now in the design phase.
Quasi-linear: In plasma theory, an expansion in powers of wave amplitude in which the fast
motion is treated linearly and zero-frequency quadratic terms determine the slow background
evolution.
Resonance: 1) A singularity in the electromagnetic plasma wave equations arising from thevariation of physical parameters in space. 2) A singularity in the particle distribution function
arising from the matching of wave and particle velocities in particular regions of phase space, awave-particle resonance.
Reynolds Stress: The advection or transport of momentum by turbulence.
RF: Radio Frequency – refers to methods of heating, flow or current drive depending on theinteraction of plasmas with externally launched waves.
ρ*: Pronounced “rho-star”, the ratio of the ion gyroradius to the minor radius. Usually ρ* ismuch smaller than unity and is used as an expansion parameter. See Gyrokinetic.
RWM: Resistive wall mode – an MHD mode driven unstable by finite wall resistivity (asopposed to zero resistivity of the ideal MHD theory).
Sawtooth: Sawtooth oscillation, a minor disruption in which core field lines reconnect flattening
the temperature profile, followed by a slow reheat of the central plasma. This occurs cyclically.
Self Organized Criticality (SOC): A characteristic of turbulent or chaotic systems where large
scale or long lasting structures arise spontaneously and are responsible for significant amounts oftransport.
Separatrix: The boundary surface between regions of a plasma on closed and open field lines.
Since there is no confinement in the direction parallel to magnetic field lines, the plasma on
open field lines is very cold (typically less than a million degrees K).
SOC: See Self Organized Criticality.
SOL: Scrape-Off Layer – region of plasma existing on open field lines outside separatrix.
Stellarator: Unlike most magnetic confinement devices, stellarators employ helical magnetic
fields formed by external coils and is intrinsically three-dimensional.
Tearing modes: Resistive MHD modes driven by plasma currents.
TFTR: Tokamak Fusion Test Reactor – a large (former) magnetic confinement experiment.

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
57

Tokamak: A magnetic confinement device relying on a large toroidal current carried by theplasma itself. (The acronym is from the Russian phrase TOroidalnaya KAmera i MAgnitnaya
Katushka, meaning toroidal chamber and magnetic coil.)
Toroidal: 1) Characteristic of a machine that is shaped like a torus or doughnut. 2)The long way
around in a torus-denoted by φ (See Fig. III.1 of the main report).
Transport Barrier: Region of strongly reduced transport - typically a bifurcated state. (See ITB,
H-mode, and Pedestal).
Transport flux: The total energy or particle flow per unit area passing through a surface
element.
Unfavorable curvature: See Curvature.
Vlasov equation: The collisionless version of the Boltzmann equation.
Wave-particle resonance: See Resonance 2).
WDM – Whole Device Modeling
Whistler wave: A high frequency plasma wave with a nonlinear dependence of frequency andwave number. Occurring naturally in the ionosphere, whistler waves were first detected by radio
operators during World War I.
X-Point: A point in the plasma cross section where the poloidal magnetic field is zero; occurs at
one or more points along the separatrix.
Zonal flows: Flows which are constant on flux surfaces but with strong radial variation. Theseflows are self-generated by plasma turbulence for which it is also an important non-linear
saturation mechanism.

FESAC ISOFS Subcommittee Final Report — Appendix
Overview of Fusion Science for the FSP
58

FRONT COVER ATTRIBUTION LIST(APPEARS ON BACK COVER)

B. Carreras, W. Dorland, L.Garcia
G. Hammett, G. Kerbel, C. Kim
J-N. Leboeuf, V. Lynch, W. ParkS. Parker, B. Rogers, C. Sovinec
D. Spong, R. Sydora, NCSX Team
W7-X Team