6200 kB PostScript reprint
N. Brummell and J. Toomre
Joint Institute for Laboratory Astrophysics, Department of
Astrophysical, Planetary and Atmospheric Sciences, University of Colorado,
Boulder, CO 80309-0440
An interdisciplinary team of researchers at several institutions (Universities of Colorado, Chicago, and Minnesota, Michigan State University, and National Center for Atmospheric Research) is working jointly on problems in geophysical and astrophysical fluid dynamics (GAFD) to utilize massively-parallel architectures to increase the spatial resolution in three-dimensional simulations of GAFD turbulence. Within our ``grand challenge'' team efforts, we are employing variously pseudo-spectral, finite-difference, multi-grid, and piecewise-parabolic method (PPM) approaches in studying the intense turbulence encountered in both planetary and stellar settings. The scale of these simulations requires corresponding progress in the computational sciences, both in order to develop and optimize software for massively-parallel computers and to capture and visualize the resulting massive data sets.
Our research related to astrophysics has so far concentrated largely on turbulent compressible convection within stars. In particular, the vigorous turbulence that results from convective instability within rotating stars not only serves to transport heat but also to redistribute angular momentum and chemical species, and can yield magnetic dynamo action. A hallmark of such turbulence constrained by rotation and stratification is that large-scale coherent structures and strong mean flows can coexist with the intense smaller-scale turbulence. Our work with three-dimensional direct simulations of turbulent compressible convection influenced by rotation is motivated by trying to understand how stars and giant planets rotate differentially, and in particular to examine the redistribution of angular momentum that leads to differential rotation within a star like the sun. As one example, these studies seek to resolve a major challenge raised recently by helioseismology, wherein acoustic oscillations of the sun are used to probe the structure and dynamics of the solar convection zone. Preliminary results from helioseismology (cf. Gough & Toomre 1991) suggest that the sun is rotating differentially with depth and latitude in a manner very different from predictions based on earlier numerical simulations. Those spherical shell simulations of rotating anelastic convection (e.g., Gilman & Miller 1986; Glatzmaier 1987) dealt with effectively laminar convective flows, whereas our recent work suggests that fully turbulent compressible flows appear to yield very different transport properties for angular momentum. This holds out the hope that in due course simulations of GAFD turbulence carried out in appropriate spherical geometries and with the inclusion of adequate stellar physics (realistic equations of state and opacities, ionization effects) may help to resolve the basic dynamical puzzle now being raised by helioseismology. However, understanding such nonlinear dynamics at a fundamental level raises formidable challenges because of the broad range of physical scales that must be resolved. High-performance computing offers the opportunity to make substantial inroads for studying the properties of such astrophysical turbulence.
The basic difficulty with trying to model turbulence in most natural
settings, and certainly so in a stellar convection zone, is that the active
dynamical scales of convection encompass a formidable range.
Turning to the sun, we expect that the largest-scale turbulent flows deep
within the convection zone should possess scales comparable to the overall
depth of that zone (occupying the outer 30 percent by radius, or about
200,000km in depth). Yet the density scale height near the top of the
zone is only about 100km, and one readily observes there intense turbulent
convection in the form of granulation with horizontal scales of order 1,000
km. Detectable there, too, are supergranulation flow patterns with scales of
order 30,000km. However, we can estimate that viscous dissipation is
operating at scales of order 0.1km or smaller, suggesting that the active
dynamical scales of turbulence in such a zone range at least from 
to
km, or thus over six decades in each physical dimension. The
most intricate three-dimensional turbulence simulations to date have just
begun dealing with order 
spatial modes, capturing three decades in
each dimension. As we will discuss below, such 
spatial zone
calculations to study the evolution of turbulence for typically 
time
steps make rather interesting demands on processors, memory, communications,
and data capture. Though we may have to wait for quite some time before
zone simulations become tractable (if ever), we can take some
comfort in realizing that the hints of structured turbulence first realized
with 
simulations are now seen in fascinating detail with the
``hero'' calculations of the day. Another three decades of
resolution (and degrees of freedom) in each physical dimension may
well herald other surprises in the resulting turbulent dynamics. The
optimism of theoretical physics encourages great extrapolations to be made,
and in that spirit we hope that the turbulence calculations of
-class likely by the end of this decade and century may provide
substantial clues about the inner workings of a real star like the sun.
Yet we recognize that all such modelling is of the large-eddy simulation
(LES) variety, requiring some manner of sub-grid-scale (SGS) representation
of the unresolved scales.
The nature of our computational challenges are revealed by considering some
of the attributes of a 
-mode turbulence simulation. For fully
compressible convection, we typically use as working variables the three
components of velocity 
, pressure p, and temperature T.
This implies that the memory requirements are 160 GB, based on 
spatial zones 
5 variables 
factor 4 of working space, for
a total of about 20 Gwords. On a 1024-node massively parallel processor
(mpp) system, we effectively need 160MB of memory per node. Turning to
estimates of operation counts for these simulations, the floating point
operation (flop) count per zone per time step range from about 200 for
pseudo-spectral codes to about 4000 for PPM codes. Since the intricate
dynamics needs to be evolved for about 100,000 times steps, the total
operation count for a given simulation are 
flop using
pseudo-spectral codes and 
flop using PPM codes. Once a
real Tflops
machine (1000 Gflops
) is realized, such a pseudo-spectral run
could be accomplished in 
s (about 6 hours), whereas the PPM
run requires 
s (about 110 hours). At the sustained 20
Gflops
that is about the best performance now, these numbers translate into
s (280 hours) and 
s (5500 hours).
The choice of computational fluid dynamics algorithms to be favored on a given mpp architecture is considerably influenced by the effective internal communication speed between processors. The use of fully explicit and spatially local methods, such as PPM, places the least demands on such communication. With PPM, our 3--D simulations on regular grids can be subdivided into many subvolumes, with only values of variables on boundary zones a few layers deep having to be exchanged after each time step between processors working on adjacent subvolumes. If the operation count is high, this is a relatively infrequent event. In contrast, pseudo-spectral codes assign successive physical planes to the processors, and carry out fast transforms on these planes within processors, followed by data transposes when calculations are required in the final dimension. Alternatively, the transforms may be carried out in parallel across processors, with either approach requiring fairly intensive data movement. Since the operation counts with pseudo-spectral codes are typically far more modest than with PPM, communication currently places major constraints upon the actual performance realized with the principal mpp contenders.
Of greater concern is how to extract and capture sufficient data to be used
for subsequent visualization and analysis of these intricate turbulent
flows. A single data snapshot of five fields and two tracers at
-resolution yields 56GB per dump. If we sample the time evolution
fairly sparsely after every 100 time steps (or 
samples per run),
then the full data set is about 56TB per run. Not only is that number
rather large given most current mass storage facilities, but the data flow
rate will also be problematic as the computing machines get faster. On
current 20Gflops
machines, the required data outflow for the PPM code is a
tractable 2.8MBs
, and for the pseudo-spectral code (with about 20-fold
fewer operations per time step) a more difficult 56MBs
. Present multi-TB
mass stores have sustained transfer rates of 3 to 20MBs
. Data output
from a Tflops
machine would in turn be 0.14GBs
(for PPM) and 2.8GBs
(for pseudo-spectral), or thus at worst a factor 
greater than
current mass store rates. Substantial improvements in rates of data output
and capture are anticipated, but it is likely that such data movement
places the severest restraints on efficient simulations of turbulence. The
options to help survive this looming data crunch are to sample full dumps
less often, or to take compressed dumps, or used weighted averages of the
data, or to do smaller problems. More acceptable is to recognize that
investments must be made for data movement and capture that may turn out to
be comparable in cost to that of the mpp itself.
Figure: Perspective full domain view of vertical velocity at one instant
in time in 3-D turbulent convection simulation (
). Here
, with the rotating f-plane
domain at latitude 
. The rendering has upward-flowing, hot
fluid as bright, with downward-flowing, cold fluid as darker tones. Note
the curvaceous network at the upper surface and the smaller-scale
turbulence at the bottom.
Original PostScript figure (1496 kB)
Our goals have been to first study the dynamics of increasingly turbulent
states in simplified ``local area'' geometries with the influence of rotation
incorporated via f--plane modelling, and then to proceed toward a more
physical ``global'' geometry encompassing an entire rotating spherical shell
of fluid. Our local-area models consist of a planar layer of fluid
positioned tangentially to a rotating sphere at a latitude 
. The
model configuration for the fully compressible treatments is essentially
that of Cattaneo et al. (1991). The physics of the model is simplified by
employing a perfect gas, and fluid properties such as the specific heats,
shear viscosity, and thermal conductivity are also assumed constant. The
flow is confined between two horizontal, impenetrable, stress-free
boundaries, and the fields are assumed to be periodic in the horizontal.
The upper boundary is held at a given temperature, whereas the heat flux is
prescribed on the lower boundary. We then proceed to solve as an initial
value problem the partial differential equations for the conservation of
mass, momentum, and energy, along with the equation of state. The
simulations are characterized by several nondimensional parameters: the
Prandtl number 
, the ratio of viscosity to thermal conductivity, the
Taylor number 
, the ratio of Coriolis to diffusive terms, the Rayleigh
number 
, the ratio of buoyancy forcing to viscous and thermal
dissipation, and the Rossby number 
, the ratio of nonlinear advective
to Coriolis forces.
Figure: Perspective 3--D view of enstrophy (vorticity squared) at one
instant in time as a companion to Figure 1. High enstrophy is bright and
opaque whereas low enstrophy is dark and translucent. Distinctive
enstrophy structures are present, with evolving cellular flows near the top
of the layer yielding sheets and tubes of vorticity, whereas in the intense
turbulence at greater depths the vorticity is often organized into
randomly aligned tubes. Large coherent tube-like structures intermingle
with small-scale random vortex tubing in the interior.
Original PostScript figure (1136 kB)
Our main computations shown here involve the solution of the fully
compressible equations using a hybrid finite-difference and pseudo-spectral
code. The vertical structure is treated by fourth-order finite differences
in the interior. The horizontal components are treated by a pseudo-spectral
method that calculates all linear operations and derivatives in spectral
space 
and performs the nonlinear multiplications in
configuration space 
, with the transform between spaces achieved
by fast transform methods. The time discretization is based on an explicit
two-level Adams-Bashforth scheme, and the thermal conduction terms are
treated implicitly with a Crank-Nicholson method.
We have studied the basic features of fully-developed rotating compressible
turbulence, with special regard to the types of flow structures,
topologies, mean-flow generation, and energy balances that may be achieved.
The simulations typically involve a density contrast across the layer of
, spanning roughly 5 pressure scale heights. Other parameters
fixed are the ratio of the specific heats, 
, and the aspect
ratio of the box, at 
. A series of runs has been made sparsely
populating the four-dimensional parameter space (
) of
the model but revealing a rich set of results. We will briefly sample some
of the characteristics of such turbulent convection.
Let us first comment on the effects of rotation on the structure and topology of the turbulent convective flows. We find that convective flows in the upper thermal boundary layer consist of a nearly laminar, cellular network, atop a fully turbulent interior, punctuated by vertically-coherent structures emanating from the upper surface (e.g., Brummell, Hurlburt, & Toomre 1993, 1995a,b). There is pronounced asymmetry between the strong, narrow downdrafts and the weak, broad upflows. This is evident in Figure 1, which shows the local-area computational domain with the vertical velocity on the visible surfaces shaded according to its strength (see also Figure 5). A plane close to the lower surface of the domain has been removed to show clearly the difference in scales between the fully-developed turbulence there and the laminar upper surface. The smoothness of the upflow regions near the surface is achieved by rapid horizontal expansion of the ascending fluid elements which will tend to smooth out any small-scale fluctuations. The cell packing near the top is irregular, and cells are born and obliterated in an intricate time-dependent fashion. The Coriolis forces of rotation have changed the characteristics of both the surface network and the turbulent interior from that of earlier nonrotating models (Cattaneo et al. 1991). The surface network has become more curvaceous, less connected and more time-dependent with the addition of rotational constraints. These changes are due to inertial oscillations of surface flows induced by the underlying rotation which mobilizes the network, and to the enhancement of a mechanism for cell destruction and creation by the background pool of vorticity. In the latter, vorticity concentrated at the interstices of the downflow network lanes is amplified by the rotational influence to produce vortices which are strong enough to evacuate their interiors, thus causing a reversal of the vertical flow and the eventual destabilization of the vortex itself. These ``spinners'' break apart the network at the junctions and allow new upflowing cells to form there.
Figure: Phase diagram of pointwise convective flux 
as function of
the vertical velocity in a horizontal plane near midlayer for turbulent
convection with 
, 
and 
, with
. The boxes A, B and C delineate horizontal sites of upflow,
downflow, and strong downflow within the horizontal sampling plane.
Original PostScript figure (1670 kB)
We find that the nature of the flows in the turbulent interior are subtly altered by rotation. The Coriolis forces provide a linear mechanism for coupling the vertical and horizontal motions, and thus with rotation the isotropy of energy and enstrophy is enhanced. This manifests itself as the concentration of vorticity into distinctive tube-like structures with a very random orientation in the interior of the flow. These vortex tubes are illustrated in Figure 2, which displays the enstrophy (vorticity squared) field in a volume rendering. Near the surface, the laminar cellular flow is dominated by vortex tubes where the coherent downflows form, and by weaker sheets of vorticity generated by the downflow lanes. At greater depths, the small-scale turbulence is characterized by intricately interwoven small intense vortex tubes. Larger structures with vertical coherence are always in existence, and their presence in conjunction with Coriolis effects can lead to enhanced Reynolds stresses and thus a generation of mean flows (cf. Brummell, Xie, Toomre, & Baillie 1995). The vortex interactions of the small-scale turbulent tubes and the larger-scale vertically-coherent strong downflowing structures provide a horizontal mixing mechanism which retards the vertical mixing by the convection, yielding a mean thermal stratification in the vertical which is substantially nonadiabatic over much of the interior. This has the important and surprising consequence that with rotation, buoyancy driving of such turbulence will be experienced throughout the layer depth, and not just dominantly in the thermal boundary layers near the top and bottom of the layer.
Figure: Horizontally-averaged convective flux 
as function of depth
for the case sampled in Figure 3, revealing that the near cancellation of
enthalpy and kinetic fluxes within the downflows leads to the upflows
dominantly carrying the convective flux.
Original PostScript figure (25 kB)
Figure: Vertical velocity w (left) and convective flux 
(right)
shown in three horizontal slices near the upper surface, middle and lower
surface in single snapshot of vigorous turbulent convection for 
, 
, 
, and 
= 2.24. Upflows are
lighter tones, downflows are darker tones. The change from smooth, nearly
laminar flows at the top to disordered, turbulent flows deeper down is
clearly evident. For 
, regions of upflow are blanked out with neutral
grey, and downflow regions encoded either as white or black according to
the sense of 
. Figure 4 indicates that horizontal averaging of these
spatially intermittent fluxes in downflows leads to near cancellation.
Original PostScript figure (1908 kB)
There is a distinctive transition from laminar to turbulent flows at
greater depths when the Prandtl number 
in the simulations was
decreased below a value of about one-third. The flows are vigorously time
dependent and turbulent, yet they exhibit coherent flow patterns
recognizable over several scale heights, contrary to the assumptions of
mixing-length approaches. Further, the surface appearance of the
convection, always being cellular and irregular, does not belie whether the
deeper structures are laminar or turbulent. The transition to turbulence
also results in remarkable changes in how the convection transports the
energy. Namely, in the laminar cases the strong downflows serve to both
carry most of the enthalpy flux upward and a significant kinetic energy
flux downward: the upflow sites make only a modest contribution to both
fluxes. However, once there is transition to turbulence at the greater
depths, there is a prominent change in the role of the downflows, in which
the (upward) enthalpy and (downward) kinetic flux contributions now nearly
cancel each other, leaving the disorganized turbulent motions outside of
those coherent structures to effect most of the overall energy transport.
Thus whereas the larger scales control the convective patterns and their
evolution, the transport is left to the smaller turbulent scales. This is
a striking example of how transports achieved by convection can be
radically modified by transitions to turbulence. Further, there are major
fluctuations in the total convective flux (sum of enthalpy and kinetic
fluxes) from site to site across any given horizontal plane, as the phase
diagram of vertical velocity w and convective flux 
in Figure 3
reveals. In the upflows (A), the pointwise 
values are sufficiently
skewed that there is a net upward flux (negative) when averaged over all
points, whereas the distribution is more nearly symmetrical for the
moderate downflows (B) and so too for the strong downflows (C), with only a
small value for the net transport emerging after averaging. This is indeed
seen in Figure 4 of the horizontally-averaged 
with depth as
contributed in turn by the upflows and the moderate and strong downflow
regions. The contributions from the upflows dominate in the convective
transport.
The results shown in Figure 4 are extraordinary, implying that throughout
the turbulent layer nearly all the convective flux is carried by the
upflowing material and not the downflows, despite the fact that the
downflows are much more coherent in space and time compared to the random
fluctuations of the smaller scale upflows. The near-zero horizontal
average value in the downflow at any depth is due to a cancellation of the
two components of the convective flux, i.e., the enthalpy flux and the
kinetic flux. Cattaneo et al. (1991) proposed a theory as to why that
might be, suggesting that cancellation of the constituent fluxes at each
point in the downflows is due to a Bernoulli-like balance achieved in the
dynamics, provided that the motion is nearly steady, the mean
stratification is nearly adiabatic, and viscous and thermal conduction are
negligible. These conditions are reasonably well met in the turbulent
simulations. Figure 5, however, points out that this theory may be overly
simplistic. Shown is the vertical velocity and measures of the convective
flux 
for three depths at a typical time in the simulation. The
convective flux is shown only in the downflows (the upflows are obscured
with a grey mask) with black marking regions of positive convective flux
and white marking regions of negative flux. This rendering reveals that
the cancellation is neither pointwise nor even local. If it conformed to
the theory, then the whole downflow regions would be colored neutrally
gray, or at least would have a very high frequency speckling between small
black and white values. Instead, we see interleaved patches of black and
white where the values of 
are actually quite large and definitely
non-zero. Surprisingly, however, the horizontal average (i.e., the sum of
all the black and white patches) is nearly zero. The cancellation of the
convective flux appears to be more a phenomenon related to major
intermittency in the flux production, both in time and space, rather than
by pointwise cancellation as was first thought. These are rather subtle
matters, emphasizing that highly turbulent flows can possess transport
properties very different from those of laminar flows.
This work was supported in part by NSF through grant ECS--9217394 and by NASA through grants NAG5--2218 and NAG5--2256. The simulations were carried out at the Pittsburgh Supercomputer Center under grant MCA93S005P.
Brummell, N. H., Hurlburt, N. E., & Toomre, J., 1995a,b, ApJ, submitted
Brummell, N. H., Xie, X., Toomre, J., & Baillie, C. 1995, in Proc. GONG'94, Helio- and Asteroseismology, ASP Conf. Ser., ed. R. K. Ulrich (San Francisco, ASP), in press
Cattaneo, F., Brummell, N. H., Toomre, J., Malagoli, A., & Hurlburt, N. E. 1991, ApJ, 370, 282
Gilman, P. A., & Miller, J. 1986, ApJ, 61, 585
Glatzmaier, G. A. 1987, in The Internal Solar Angular Velocity, eds. B. R. Durney & S. Sofia (Dordrecht, Reidel), p. 263
Gough, D. O., & Toomre, J. 1991, ARA&A, 29, 627
6200 kB PostScript reprint