FAQ on "Flux-driven" vs. "Gradient-driven" turbulence

             Greg Hammett, PPPL, March 7, 2002

From time to time I am asked questions about "flux-driven"
vs. "gradient-driven" turbulence simulations.  Sometimes it is phrased
in terms of "constant-flux" boundary conditions vs. "periodic" boundary
conditions.  Below is a lengthy email I wrote to a colleague on this
topic.

------------------------------

The following text is rather long.  To try to put the issue succinctly:


Occastionally I am asked questions about "flux-driven"
vs. "gradient-driven" turbulence simulations.  Sometimes the questions
are phrased in terms of "constant-flux" boundary conditions
vs. "periodic" boundary conditions.  Some people think that our
flux-tube simulations with periodic blundary conditions have forced the
gradient to be artificially fixed and that we therefore can't see
avalanches or similar effects due to fluctuating gradients.  They
somehow seem to think that our temperature gradient is everywhere a
constant, or at the edges is constant, or something like that.  But
these are misconceptions.  At any radial position in our simulations,
the temperature gradient is allowed to self-consistently fluctuate in
response to the turbulence.  It is only the temperature gradient
averaged over the whole box which is assumed constant in a periodic
simulation, and for a large enough simulation this assumption should be
fine.  Not only can we see avalanches, we can let them propagate a very
long distance if they want to, propagating out one side of the box and
entering back through the other side.  

Both "constant-flux" and "periodic" boundary condition codes should
converge to the same answer for sufficiently large box size (if properly
implemented).  For reasons described below, I believe periodic boundary
conditions are more efficient and will converge more quickly (for
homogeneous systems where they are appropriate).

------------------------------

Here are some answers to the questions you asked.  This is much longer
than you (or I) thought it would be, and you already know much of this.
But I've been asked similar questions several times by various people,
and several misconceptions seem to circulate, so I wanted to write down
an answer that I could refer a wide range of people to.  I welcome any
feedback or further comments from you.

Before I answer your questions about "periodic" vs. "constant-flux"
boundary conditions, let me first summarize the standard motivation for
using periodic boundary conditions.

Periodic boundary conditions are natural boundary conditions when
dealing with small-scale homogeneous turbulence embedded in a larger
system.  Ron Waltz has articulated the standard motivation for periodic
boundary conditions by the following logic. If a simulation domain is
large enough, then, for a statistically homogeneous system (independent
of position), the fluctuations at one position should statistically look
the same as the fluctuations at another position.  I.e., the range of
eddy shapes and sizes seen exiting or entering one side of the box
should be the same as the range of eddy shapes and size seen entering or
exiting the other side of the box.  If these two opposite sides of the
box are sufficiently far apart so that they are uncorrelated, we can
satisfy this statistical equivalence by imposing exact equivalence.  In
other words, for every eddy of a certain size and shape which exits the
box on one side, exactly the same eddy enters the other side.  This is
what the assumption of periodicity accomplishes.  If the two sides of
the box are sufficiently far apart, this doesn't matter.

As you know, periodic boundary conditions have long been used for
small-scale homogeneous turbulence simulations (such as Orszag and many
others do), in simulating Rossby wave turbulence in planetary flows, in
Hasegawa-Mima/Terry-Horton-Waltz types of simulations, etc.

Interestingly, there has been a recent discovery of a very important
type of plasma instability that drives turbulence in accretion disks.
The first nonlinear simulations of this problem were with a "shearing
box" geometry, which is very similar to the tricks we do in flux-tube
simulations to properly impose radial periodicity in the presence of
magnetic shear, or that Dimits and Waltz used to handle the effects of
sheared flows in periodic flux-tube simulations.  See for example:

"Instability, Turbulence, and enhanced transport in accretion disks",
S.A.  Balbus, J.F. Hawley, Rev. of Modern Physics 70, 1 (1998).

Of course there are times when periodic boundary conditions aren't
appropriate: if the turbulence is at sufficiently large scales that
interactions with the system-size and system-dependent geometry are
important (i.e., the homogeneous assumption breaks down), or if there
are other sources of significant inhomogeneities in some direction.  In
fact, it's interesting that periodic b.c.'s are used to study the
inertial range of homogeneous Kolmogorov turbulence, even though the
Kolmogorov spectrum E(k) ~ k^(-5/3) is always peaked at long wavelengths
(or at least at the stirring scale).  Technically this means that the
box-size scale (or the scale at which stirring occurs) is important, but
they are usually interested in the statistics at the turbulence at the
smaller scales of the inertial range, and the assumption is that the
large scales don't impact that much.

Attached are a few drawings (in constant-flux.pdf) to illustrate some of
the points made below.

-----------------------------------------

Now let's turn to your specific questions:

> 
> Hi Greg,
> hope you're well in London.
> I've got a little question:
> 
> ... you once mentioned to me that you think
> constant flux boundary conditions (... "flux driven transport" ...)
> do not produce any new phenomena compared to the usual periodic
> boundary conditions in a flux tube.  Do you have a convincing argument
> on this, in particular one that clearly disproves the line of thought
> below?
> 
> Somebody in favor of the flux driven transport might argue like this:
> Consider a system with a relatively hard gradient threshold for
> transport (typical example: ITG). Applying constant gradient boundary
> conditions will not produce any transport until the gradient exceeds
> the threshold, where the transport shoots up.
> 
> With constant flux boundary conditions, applying only a minor amount
> of heat flux (by sufficiently spatially smooth sources within the
> computational domain, which is supposed to have impenetrable
> boundaries) will raise the gradient, until its threshold, where the
> instabilities will begin to grow. This growing will take some time, so
> that the gradient will be able to rise above the critical value. At
> this point the transport sets in pretty strongly, pushing the gradient
> way below the critical value.

If the box is wide enough radially, it is very unlikely to push the
gradient below the critical value everywhere at the same time.  Instead,
at some radii the temperature gradient will be above the critical value
(temporarily), and at other radii the temperature gradient will be below
the critical gradient.  The temperature gradient dT(r,t)/dr at any
position will fluctuate in time, sometimes being above the critical
gradient and sometimes being below the critical gradient (maybe even
sometimes reversing sign).  If a region of strong temperature gradient
(above the critical gradient) is seen to propagate radially significant
distances, one might term that an "avalanche".  By an ergodic theorem
that holds for more systems of interest, the time-average of the
temperature gradient will be the same as the spatial-average of the
temperature gradient (if the box is large enough), so that the
temperature gradient in the constant-flux code will fluctuate around the
same average temperature gradient used in a periodic code for the same
level of average heat flux.

Simulations with either "radially periodic" boundary conditions or with
"constant flux" boundary conditions (if properly constructed) are both
be able to see this behavior of fluctuations in the local temperature
gradients and local heat fluxes.  Both types of boundary conditions
should give the same plot time-and-space averaged heat flux
vs. time-and-space averaged temperature gradient.

There is a misconception that some people have that a code with periodic
boundary conditions can't see behavior such as avalanches because they
"fix the temperature gradient" instead of "fixing the flux".  But this
simply isn't true: only the average temperature gradient is assumed
fixed while at any given radial position the temperature gradient
fluctuates self-consistently as heat is turbulently moved around.
Periodic boundary conditions do indeed allow avalanches, and in fact
allow avalanches to propagate a very long distance, as they will allow
an avalanche to propagate out one side of the box and back in the other
side.  In that sense, periodic bc's are actually better than constant-flux
bc's, which prevent the propagation of anything across the boundaries
and have buffer regions near the edge where artificial dynamics must be
introduced to model what happens in adjacent regions.

Of course either approach should converge to the right answer in the
limit of a sufficiently large box, which is why in practice it is useful
to do box-size convergence tests from time to time.  (I think that
periodic boundary conditions will converge to the right answer more
quickly than constant flux boundary conditions, for reasons described
more below.)

>  The turbulence will now still go on a
> little, because it takes some time to damp out, so that after it has
> "died", the gradient is definitely subcritical. Then the game starts
> again, we will get a cyclic oscillation in the gradient, with the
> possibility that the average gradient is actually below the critical
> gradient.

There are two ways in which such cases might occur.

One way is to consider isolated avalanches, isolated regions of steep
temperature gradient separated by wide regions where the temperature
gradient is below critical.  Clearly, in order to properly handle such a
case, a code with periodic boundary conditions would need to be wide
enough to contain at least a few avalanches (regions of steep
temperature gradient).  In practice, I've never seen regions of steep
temperature gradient that are separated by large distances in anyone's
simulations, so the box doesn't need to be very wide to converge to the
right answer.  

The other way to look at this is more like you describe: perhaps regions
of super and sub-critical gradients are not far apart, but the
time-averaged temperature gradient is below the critical gradient.

Such sub-critical turbulence can in fact be studied in codes with
periodic boundary conditions, it just requires initializing the
turbulence to a sufficiently high amplitude, or starting with the
turbulence from a super-critical case and then reducing the temperature
gradient to subcritical.  One can see evidence for sub-critical
turbulence in a code with periodic b.c.'s by plotting chi vs. R/L_T
(where L_T is the average gradient scale length) for a range of R/L_T
above the linear threshold.  If sub-critical turbulence exists, then one
should find that chi vs. R/L_T extrapolates to non-zero chi for R/L_T <
R/L_Tcrit.  Bruce Scott and Jim Drake and Barrett Rogers have in fact
studied such finite-amplitude sub-critical turbulence for edge
drift-wave parameters. (Drake and Rogers did it with radially periodic
simulations.  I think Bruce might have done it with "constant flux"
buffer zones or with both boundary conditions, I don't remember, but in
any case Bruce didn't make a big deal out of the boundary conditions as
being important for the sub-critical behavior.)  But for core
ITG/ETG/trapped-electron turbulence we've never seen any indications of
subcritical turbulence.  If anything, core turbulence is supercritical
because of the Dimit's nonlinear shift (the turbulent generation of
zonal flows that suppress transport near marginal stability)!

> 
> I.e., with constant flux BC it is possible to get transport with the
> gradients below the critical value, with constant gradient it is
> not. Also, the cyclic behavior can be interpreted as some kind of
> intermittency, which is absent in the other case. It is now no quantum
> leap of thought to argue, that indeed the tokamak is fueled rather by
> constant sources, and not by constant gradients (heating neutral
> fueling).
> 

In fact, this intermittent behavior is what shows one of the
inconsistencies of the "constant flux" boundary conditions.  You can try
to put in heat in the buffer zone at the left boundary at a constant
rate and take it out in the buffer zone at the right boundary at a
constant rate, but the flux of heat across the middle of the box will
fluctuate in time, at least some.  This is what the plasma wants to do,
and the simulations with periodic b.c.'s allows the plasma do this
self-consistently if it wants to.

It is true that the external heating source is fairly uniform in time,
but that is primarily concentrated near the core region and it does not
follow that the heat flux is uniform in time anywhere.  Further out, the
plasma is heated by heat that turbulently flows from deeper inside the
plasma.  Even in the core region, the heat flux may fluctuate, since the
energy balance equation is

(3/2) n dT/dt = P(r) - Div(Q)

and even if the heating power P(r) is independent of time, the heat flux
across any surface may fluctuate as dT/dt fluctuates.

In short, I don't see any advantage to a "constant flux" boundary
condition over a periodic boundary condition for small scale turbulence.
The argument that one needs a constant flux near the boundaries because
temperature gradient might be fluctuate is self-contradictory, since the
temperature gradient can only fluctuate in time if the divergence of the
flux is not constant in time.

(Of course, if the flux surface of the real plasma or of the simulations
is sufficiently large compared to the eddy sizes, then averaging over
the surface will average over many eddies and will result in a flux
which fluctuates less in time.  A subtlety is that zonal flow dynamics
bring in the dynamics of the n=m=0 mode, which has been found to make
the flux across any given surface fluctuate a fair amount even after
averaging over the whole flux surface.  Klaus Hallatschek (and Ben
Carreras pointed out something similar) has pointed out that, in
principal, if the box is made large enough, then the role of a single
mode (the n=m=0 zonal mode) can't matter anymore.  Part of the reason
the simulations are fairly insensitive to box size, even for the
relatively small sizes typically used in flux-tube simulations, and even
when zonal modes are important, is that it turns out that the zonal
modes are actually driven not just by random turbulent noise but by a
nonlinear secondary instability whose growth rate is enhanced by the
peculiar adiabatic electron response for n=m=0 modes.  This is
essentially a Kelvin-Helmholtz instability that has also been analyzed
in terms of a modulational instability (see papers by Drake, Guzdar,
et.al., Diamond et.al., Liu Chen et.al., Rogers and Dorland et.al.).  It
is probably still true that for an extremely large system the zonal
modes will become unimportant, but that would require parameters such
that the electrons can't transit around a flux-surface quickly enough in
a turbulence-decorrelation time to set up the proper adiabatic electron
response that enhances the growth rate of the zonal modes.  For typical
ITG/trapped-electron core parameters I think the present simulations are
okay.)


> I'd be grateful for your insights.
> (I have myself the strong feeling the constant flux BCs are bogus,
> however I'd be very interested in your argument).
> 

Both "constant-flux" bc's and "radially periodic" bc's should converge
to the same answer for a large enough box.  Constant-flux boundary
condition codes have to implement some kind of buffer region near the
radial edges of the simulation to provide sources and sinks and
artificially damp out the potential.  (This is part of why I think
"constant-"flux" boundary conditions is something of a misnomer, as
there is no simple way of fixing the flux across a boundary and they
have to alter the equations in a buffer zone of at least several
gyroradii or banana widths near the boundaries.)  They have to calculate
the statistics of the turbulence by averaging over a volume which is
more than a radial decorrelation length away from these artificial
buffers, which wastes even more simulation space.  Thus I think a code
with radially periodic boundary conditions will converge faster because
it doesn't need to waste space on these artificial buffer regions.

If one looks at the flux-surface-averaged temperature profile T(r,t), in a
"constant-flux BC" code, it can always be written in the following form:

T(r,t) = T_0(t) ( 1 - r/L_T(t)) + sum_k T_k(t) exp(i k r)

I.e., any arbitrary function over some domain can be written in terms of
a linear piece and periodic pieces.  The only difference between a
periodic and a constant-flux code is that the spatially-averaged
temperature T_0(t) and the spatially-averaged temperature gradient scale
length L_T(t) are technically allowed to fluctuate some in a
constant-flux b.c. code, while they are assumed constant in a periodic
b.c. code.  If the box is sufficiently wide radially, then in fact T_0
and L_T will be relatively constant, thus justifying the periodic
b.c. assumption.

Nevertheless, there are cases where the periodic b.c. assumption breaks
down, such as if the radial non-homogeneity of the plasma is becoming
important.  I.e., if the simulation is wide enough that omega_*(r) or
eta(r) or other parameters vary enough radially so that the
characteristics of the turbulence at r_1 and r_2 (the two edges of the
simulation) are not statistically identical.  Certain types of
inhomogeneity can still be modeled with a suitable modification of the
periodicity constraints.  Even the shear in the magnetic field required
a modification of simple radial periodicity, which led to the twisting
periodicity used by Beer, Cowley, and Hammett, by Dimits, by Scott, etc.
Sheared equilibrium-scale flows can also be included in
modified-periodic simulations by effectively using a moving sheared box
as done by Dimits et.al. or by Waltz et.al.  (this is like the
Balbus-Hawley shearing-box simulations used for astrophysical accretion
disk simulations).  But for stronger radial inhomogeneities that can't
be handled by twists or shifts in the periodic conditions, one is then
forced to some kind of radial boundary condition that doesn't require
periodicity.  This can be particularly important in relatively small
tokamaks (but becomes less important at reactor scales) where shear in
omega_*(r) can have a stabilizing influence similar to shear in the
equilibrium-scale poloidal ExB velocity.  Waltz's APS 01 invited talk
discusses some of their work on this employing non-periodic b.c.'s,
though I believe his conclusion is that local flux-tube simulations with
periodic b.c.'s give similar results (unless the shear omega_*(r)
etc. is extremely strong), once he included some model source/sinks in
the non-periodic simulation to prevent artificial relaxation of the
long-time temperature profile.  Non-periodic radial boundary conditions
are of course useful for edge turbulence simulations that explicitly
include open and closed field lines, losses to the divertor plate,
and/or strong radial variation in omega_*(r), etc., where radial
inhomogeneities are important, etc.

A side note:

While 1-D sand-pile models and the "avalanche" paradigm are interesting
to study, in practice I don't think they are very relevant to
small-scale 3-D fusion plasma turbulence, at least away from the edge.
I describe some things above in terms of avalanches just because that is
what some people seem to have in mind, and to emphasize that periodic
boundary conditions do not exclude them.  Or perhaps there are
avalanche-like events at scales of order the turbulent decorrelation
length, but we are interested in statistics averaged over larger scales
than that.  Even if there are cases where the sand-pile avalanche
paradigm is relevant, then they would still need to be studied in direct
nonlinear 3-D simulations of primitive gyrokinetic/gyrofluid equations
to determine the parameters that would go into a simpler avalanche
model.

The classic simple sand-pile model has avalanche events of sizes (with a
"1/f" distribution in time and space) so one way to look for evidence of
avalanches is simply to increase the box size and see if the transport
continues to go up.  Since we usually see that the transport eventually
becomes independent of box size for sufficiently large box size, I
conclude that there is some long-wavelength cutoff to the size of
relevant eddies or avalanches.  (Of course, near the edge of a tokamak
the turbulence scale lengths become comparable to the gradient scale
lengths and everything is non-local and harder to treat accurately.)  If
the transport continued to increase as the box size increased, one might
argue that this was evidence for avalanche or streamer behavior, so a
local flux-tube simulation is breaking down.  But one should note that
"constant-flux" boundary conditions don't help in this situation, since
one would also find the transport increasing with box size in that type
of simulation if real 1/f avalanches existed.  In such a case, you would
really have to simulate the full device size, resolving all scales
simultaneously.  T.S. Hahm's 2001 APS invited talk found no evidence of
non-diffusive transport in Z. Lin's global simulations.  Also,
experimentally measured temperature profiles are fairly constant in time
on large scales, and have only small amplitude fluctuations on
relatively small time and spatial scales (except for things like
disruptions or MHD instabilities which we know have their origin in
other phenomena).