A 15.2 TFlops Simulation of Geodynamo on the Earth Simulator
Akira Kageyama,∗ Masanori Kameyama, Satoru Fujihara, Masaki Yoshida, Mamoru Hyodo, and Yoshinori Tsuda
Earth Simulator Center, Japan Agency for Marine-Earth Science and Technology, Yokohama 236-0001, JapanFor realistic geodynamo simulations, one must solve the magnetohydrodynamic equations to follow time development of thermal convection motion of electrically conducting ﬂuid in a rotating
spherical shell. We have developed a new geodynamo simulation code by combining the ﬁnite difference method with the recently proposed spherical overset grid called Yin-Yang grid. We achieved
performance of 15.2 Tﬂops (46% of theoretical peak performance) on 4096 processors of the Earth
Simulator.

I.

INTRODUCTION

The magnetic compass points to the north since the Earth is surrounded by a dipolar magnetic ﬁeld. In the 19th
century, Gauss showed that the origin of this geomagnetic ﬁeld resides in the Earth’s interior. It is now broadly
accepted that the geomagnetic ﬁeld is generated by a self-excited electric current in the Earth’s core, which is adouble-layered spherical region of radius r = 3500km. The inner part of the core, called the inner core (radius
1200km), is iron in solid state, and the outer part, called the outer core (radius r: 1200km ≤ r ≤ 3500km), is also
iron but in liquid state due to the high temperature of the planetary interior. The electrical current is generated
by magnetohydrodynamic (MHD) dynamo action—the energy conversion process from ﬂow energy into magnetic
energy—of the liquid iron in the outer core. This geodynamo phenomenon is the target of our simulation reported inthis paper.
In the last few decades, computer simulation has emerged as a central research method for geodynamo study [10].
From the beginning of this new and exciting simulation area, we have been making important contributions by
combining large scale simulation and advanced visualization technology: demonstration of the strong magnetic ﬁeld
generation by MHD dynamo [4], physical mechanism of the dipole ﬁeld generation [6], and spontaneous and repeatedreversals of the dipole moment (north-south polarity) [5, 11, 13].
Recently, we have proposed a new grid system for spherical geometry based on the overset (or Chimera) grid
methodology. We have applied this new spherical grid to geodynamo simulation. In this paper, we report thedevelopment of this new geodynamo simulation code and its performance on the Earth Simulator [2, 16, 20], of Japan
Agency for Marine-Earth Science and Technology (JAMSTEC), Japan.

II.

YIN-YANG GRID

Since the Earth’s outer core is a spherical shell region between the inner core and the mantle, we need a numericallyeﬃcient spatial discretization method for a spherical shell geometry in order to achieve high performance on massively
parallel computers.
In our previous geodynamo simulations, we used a ﬁnite diﬀerence method based on the latitude-longitude gridin spherical polar coordinates with the radius r (ri ≤ r ≤ ro ), colatitude θ (0 ≤ θ ≤ π), and the longitude φ(−π < 0 ≤ φ ≤ π). Due to the existence of the coordinate singularity and grid convergence near the poles of thelatitude-longitude grid, we had to take special care at the poles and this inevitably degraded the numerical eﬃciency
and the performance of our previous dynamo code.
Since the ﬁnite diﬀerence method is suitable for massively parallel computers, especially massively parallel vector
supercomputers like the Earth Simulator, we have decided to further exploit the possibilities of the ﬁnite diﬀerence
method in spherical geometries, and to explore a new spherical grid system as the base of the spherical ﬁnite diﬀerence
method.
There is no grid mesh that is orthogonal all over a spherical surface and, at the same time, free of coordinate
singularity or grid convergence problems. So, we decompose the spherical surface into subregions. The decomposition,or dissection, enables us to cover each subregion by a grid system that is individually orthogonal and singularity-free.

∗ Electronic address: kage@jamstec.go.jp; 0-7695-2153-3/04 $20.00 (c) 2004 IEEE

2

FIG. 1: Basic Yin-Yang grid. The Yin grid and Yang grid are combined to cover a spherical surface with partial overlap.

The dissection of the computational domain generates internal borders between the subregions. In the overset grid
approach [1], the subdomains are permitted to partially overlap one another on their borders. The overset grid is also
called an overlaid grid, or composite overlapping grid, or Chimera grid [17]. It is now one of the most important gridtechniques in computational aerodynamics.
Recently, we have proposed a new overset grid for spherical geometry [7]. The grid is named “Yin-Yang grid” after
the symbol for yin and yang in the Chinese philosophy of complementarity. The Yin-Yang grid is composed of two
identical and complemental component grids. Compared with other possible spherical overset grids, the Yin-Yang gridis simple in its geometry as well as metric tensors. A remarkable feature of this overset grid is that the two identicalcomponent grids are combined in a complemental way with a special symmetry. The Yin-Yang grid has already been
applied to a mantle convection simulation in a spherical shell geometry with detailed benchmark tests [21]. The YinYang grid has also been applied to simulations of the global circulation of the atmosphere, ocean, and their coupledsystem [3, 9, 14, 18, 19].
There are several varieties of Yin-Yang grid [7]. The most basic type, adopted in the geodynamo simulation reported
in this paper, is shown in Fig. 1. It has two components—Yin and Yang—that are geometrically identical (exactly
the same shape and size); see Fig. 1(a). We call the two component grids the “Yin grid” (or n-grid) and the “Yanggrid” (or e-grid). They are combined to cover a spherical shell with partial overlap on their borders as shown in
Fig. 1(b). Each component grid is, in fact, a part of the latitude-longitude grid on a spherical surface. (It is piled
up in the radial direction). A component grid, say the Yin grid, is deﬁned as a part of low-latitude region of theusual latitude-longitude grid: 90◦ around the equator in latitude (between 45◦ N and 45◦ S), and 270◦ in longitude.
Another component grid, the Yang grid, is deﬁned in the same way but in diﬀerent spherical coordinates. The virtual
north-south axis of the Yang grid is located on the equator of the Yin grid’s coordinates. The relation between Yin
coordinates and Yang coordinates is written in Cartesian coordinates by(xe , y e , z e ) = (−xn , z n , y n ),or (xn , y n , z n ) = (−xe , z e , y e ),(1)

where (xn , y n , z n ) is the Yin grid’s Cartesian coordinates and (xe , y e , z e ) is Yang grid’s. Note that the forward andinverse transformations between Yin and Yang are given in the same form, which is a reﬂection of the complementary
relation between the Yin and Yang grids.
It should be stressed that the Yin grid and Yang grid are identical in shape, geometry, grid, and metric tensors.
This means that all subroutines designed for, say the Yin grid, can be used for the Yang grid with no modiﬁcation.
In addition to this geometrical symmetry, the Yin grid and Yang grid are in complementary coordinates relative toeach other as indicated by eq. (1). This means that any interaction or communication from a grid point on Yin toa grid point on Yang is exactly the same as that from Yang to Yin. This fact makes the Yin-Yang-based computercode very concise and eﬃcient.

3
Another advantage of the Yin-Yang grid resides in the fact that the component grid (Yin or Yang) is nothing but
a (part of) latitude-longitude grid. This means that we can directly deal with the equations to be solved with thevector form in the usual spherical polar coordinates; components of a ﬂow vector v, for example, can be written as{vr , vθ , vφ } in the program. The analytical forms of metric tensors such as for Laplacian in spherical coordinates arewell known. Since we can directly code the basic equations as they are formulated in spherical coordinates, we canmake use of various resources of mathematical formulas, program libraries, and tools that have been developed inspherical polar coordinates.
Following the general overset methodology [1], interpolations are applied on the boundary of each component gridto set the boundary values, or internal boundary conditions. When one sees the component grid of the basic YinYang grid shown in Fig. 1 in the Mercator projection, it is a rectangle; the four corners intrude most into the othercomponent grid. Even if the grid mesh is taken to be inﬁnitesimally small, the overlapping area has still non-zero ratio
of about 6% of the whole spherical surface. We obtain two solutions—one from the Yin side and other from the Yang
side—in the overlapped area. This might be a slight (6%) waste of computational time, but the “double solution”
causes no problem in actual calculations. The diﬀerence between the two solutions is within the discretization errorthat is omnipresent on the sphere in any case. The diﬀerence between the solutions is so small that we do not needto “blend” the double solutions there in the post processing for the data visualization; we just choose one of the twosolutions and the resulting visualization shows smooth pictures. There is no indication of the internal border between
the Yin and Yang grids (e.g., Fig. 2 below).
If one still desires to minimize the overlapped area, it is possible by modifying the component grid’s shape fromthe rectangle. It is obvious that a Yin-Yang grid with a minimum overlap region can be constructed by a division, or
dissection, using a closed curve on a sphere that cuts the sphere into two identical parts. There are inﬁnite numberof such dissections of a sphere. Two examples of them—“baseball type dissection” and “cube type dissection” can befound in [7].
III.

BASIC EQUATIONS AND NUMERICAL METHODS

Since the physical model for the geodynamo simulation is the same as our previous one and described in ourpapers [e.g., 4], here we describe it only brieﬂy. We consider a spherical shell vessel bounded by two concentric
spheres. The inner sphere of radius r = ri denotes the inner core and the outer sphere of r = ro denotes the coremantle boundary. An electrically conducting ﬂuid is conﬁned in this shell region. Both the inner and outer sphericalboundaries rotate with a constant angular velocity Ω. We use a rotating frame of reference with the same angularvelocity. There is a central gravity force in the direction of the center of the spheres. The temperatures of both theinner and outer spheres are ﬁxed; hot (inner) and cold (outer). When the temperature diﬀerence is suﬃciently large,
a convection motion starts when a random temperature perturbation is imposed at the beginning of the calculation.
At the same time an inﬁnitesimally small, random “seed” of the magnetic ﬁeld is given.
The system is described by the following normalized MHD equations:∂ρ= −∇ · f ,
∂t

(2)

1
∂f
= −∇ · (vf ) − ∇p + j × B + ρg + 2ρv × Ω + µ(∇2 v + ∇(∇ · v)),
∂t
3

(3)

∂p
= −v · ∇p − γp∇ · v + (γ − 1)K∇2 T + (γ − 1)ηj2 + (γ − 1)Φ,
∂t

(4)

∂A
= −E,
∂t

(5)

with
p = ρT,2

g = −g0 /r r̂,

B = ∇ × A,
j = ∇ × B,
Φ = 2µ eij eij − 13 (∇ · v)2 ,E = −v × B + ηj,


∂vj
∂vi
+
eij = 12 ∂x
∂xi .
j

(6)

Here the mass density ρ, pressure p, mass ﬂux density f , magnetic ﬁeld’s vector potential A are the basic variablesin the simulation. Other quantities; magnetic ﬁeld B, electric current density j, and electric ﬁeld E are treated as

4subsidiary ﬁelds. The ratio of the speciﬁc heat γ, viscosity µ, thermal conductivity K and electrical resistivity η areassumed to be constant. The vector g is the gravity acceleration and r̂ is the radial unit vector; g0 is a constant. We
normalize the quantities as follows: The radius of the outer sphere ro = 1; the temperature of the outer sphere T (ro )= 1; and the mass density at the outer sphere ρ(ro ) = 1.
The spatial derivatives in the above equations are discretized by second-order central ﬁnite diﬀerences in sphericalcoordinates: (r, θ, φ). The fourth-order Runge-Kutta method is used for the temporal integration.
There are six free parameters—including three dissipation constants; viscosity µ, thermal conductivity K, andelectrical resistivity η—in this MHD system. We basically set the same parameter values adopted in our previous
simulation that showed repeated dipole reversals [11] except for the three dissipation constants; we set each of them 10times smaller than in the previous simulation. This means that the ﬂuid’s Reynolds number and magnetic Reynolds
number are both 10 times larger, the Rayleigh number—an index of the vigor of thermal convection—is 100 times
larger, 3 × 106 , and the Ekman number is 2 × 10−5 . Thus we can achieve more turbulent, and therefore more realistic
simulations compared with our previous simulations.

IV.

PARALLEL PROGRAMMING AND PERFORMANCE ON THE EARTH SIMULATOR

We have developed a new geodynamo simulation code (hereafter “yycore” code) for the Earth Simulator by converting our previous geodynamo code, which was based on the traditional latitude-longitude grid, into the Yin-Yang
grid. We use the Fortran90 language and the MPI library. We make use of the module unit of Fortran90 to keep
concise and neat program structure. All subroutines and functions are module subprograms. (There are no external
subprograms.) In order to realize data hiding and encapsulation in the module level, we set the default attribute ofmodule data and subprograms be private. The yycore code has about 6,200 lines of source (including comments).
There are 12 modules and 1 main program. The total number of module-subprograms are 153. Among them, 48
subprograms are public (global scope).
We have found that the code conversion from our previous latitude-longitude-based code into the new Yin-Yangbased code is straightforward and rather easy. This is because most of the Yin-Yang grid code shares source lines
with the latitude-longitude grid code: Our previous geodynamo code was basically a ﬁnite-diﬀerence MHD solver on
spherical coordinates with a full span of colatitude (0 ≤ θ ≤ π) and longitude (−π < φ ≤ π); on the other hand,
the Yin-Yang grid code (yycore) is also a ﬁnite-diﬀerence MHD solver on the spherical coordinates, but with just thesmaller span of colatitude (π/4 ≤ θ ≤ 3π/4) and longitude (−3π/4 ≤ φ ≤ 3π/4). The major diﬀerence is the newboundary condition (interpolation) for the communication between Yin grid and Yang grid.
Our experience with the rapid and easy conversion from latitude-longitude code into Yin-Yang code should be
encouraging for others who have already developed codes that are based on latitude-longitude grids in the sphericalcoordinates, and who are bothered by numerical problems and ineﬃciency caused by the pole singularity. We would
like to suggest that they try the Yin-Yang grid.
Since the Yin grid and Yang grid are identical, dividing the whole computational domain into a Yin grid part and a
Yang grid part is not only natural but also eﬃcient for parallel processing. In addition to this Yin-and-Yang division,
further domain decomposition within each grid is applied to achieve high performance on the Earth Simulator, amassively parallel computer.
The Earth Simulator, whose hardware speciﬁcations are summarized in Table I has three diﬀerent levels of parallelization: vector processing in each arithmetic processor (AP); shared-memory parallelization by 8 APs in eachprocessor node (PN); and distributed-memory parallelization by PNs.
In our yycore code, we apply vectorization in the radial dimension of the three-dimensional (3D) arrays for physicalvariables. The radial grid size is 255 or 511, which is just below the size (or doubled size) of the vector register of the
Earth Simulator (256) to avoid bank conﬂicts in the memory [16].
We use MPI both for the inter-node (distributed memory) parallel processing and for the intra-node (sharedmemory) parallel processing. This approach is called “ﬂat-MPI” parallelization.
As we mentioned above, we ﬁrst divide the whole computational domain into two identical parts—here we call them“panels”—that correspond to the Yin grid’s domain and the Yang grid’s domain in Fig. 1(a). In the yycore code,
MPI COMM SPLIT is invoked to split the whole processes into two groups. (The total number of processes is even.) The
world communicator is stored in a variable named gRunner%world%communicator, where gRunner is a nested derived
type.
For further parallelization within a panel, there are several choices for the method of domain decomposition. Herewe choose two-dimensional decomposition in the horizontal space, colatitude θ and longitude φ, in each panel. For
this purpose, we call MPI CART CREATE to make a two-dimensional process array with optimized rank order.
For the intra-panel communication, MPI SEND and MPI IRECV are called between nearest neighbor processes. Eachprocess has four neighbors (north, east, south, and west). The rank numbers for the neighbors, obtained by invoking

5Peak performance of arithmetic processor (AP)
8 GﬂopsNumber of AP in a processor node (PN)
8
Total number of PN
640
Total number of AP
8 AP × 640 PN = 5120
Shared memory size of PN
16 GB
Total peak performance 8 Gﬂops × 5120 AP = 40Tﬂops
Total main memory
10TB
Inter-node data transfer rate
12.3 GB/s × 2
TABLE I: Speciﬁcations of the Earth Simulator.

MPI CART SHIFT with the panel’s communicator gRunner%panel%communicator, are also stored in gRunner.
Communication between two groups (Yin and Yang) is required for the overset interpolation. This communication
is implemented by MPI SEND and MPI IRECV under gRunner%world%communicator.
The best performance of yycore code with this ﬂat MPI parallelization is 15.2 Tﬂops. This number was obtained
from the hardware counter for ﬂoating point operations on the Earth Simulator, which is automatically reported by
setting the environment variable MPIPROGINF. This performance is achieved by 4096 processors (512 nodes) with thetotal grid size of 511(radial) × 514(latitudinal) × 1538(longitudinal) × 2(Yin and Yang). Since the theoretical peak
performance of 4096 processors is 4096 × 8 Gﬂops = 32.8Tﬂops, we have achieved 46% of peak performance in thiscase. The average vector length is 251.6, and the vector operation ratio is 99%. These data can be seen in the output
of MPIPROGINF shown in List I. The high performance of the yycore code is a direct consequence of the simple andsymmetric conﬁguration design of the Yin-Yang grid: It makes it possible to minimize the communication time (10%)
between the processes in the horizontal directions, and enables optimum vector processing (with 99% of operation
ratio) in the radial direction in each process.
MPI Program Information:========================
Note: It is measured from MPI_Init till MPI_Finalize.
[U,R] specifies the Universe and the Process Rank in the Universe.
Global Data of 4096 processes:=============================RealTime (sec)UserTime (sec)
System Time (sec)
Vector Time (sec)Instruction Count
Vector Instruction Count
Vector Element CountFLOP CountMOPS
MFLOPS
Average Vector Length
Vector Operation Ratio (%)
Memory size used (MB)Min [U,R]

Max [U,R]:
452.157 [0,2623]
454.266 [0,680]
453.457:
441.499 [0,1741]
447.001 [0,7]443.220:4.232 [0,455]
5.499 [0,9]4.498:
321.969 [0,245]
380.540 [0,2138]
351.678: 45628623153 [0,3152] 48004139033 [0,2]46732455581: 13505506552 [0,855]
14116925576 [0,2178] 13758270302: 3395555533860 [0,1904] 3553256715628 [0,2178] 3461109543510: 1641123079258 [0,52]
1668645286059 [0,2178] 1642792822350:
7706.674 [0,7]8102.270 [0,3938]
7883.403:
3671.412 [0,7]3769.202 [0,1890]
3706.499:
250.582 [0,1120]
252.856 [0,4068]
251.564:99.000 [0,7]99.105 [0,3135]
99.056:1042.944 [0,0]1122.913 [0,9]1106.882

Overall Data:=============RealTime (sec)UserTime (sec)
System Time (sec)
Vector Time (sec)GOPS(rel. to User Time)
GFLOPS (rel. to User Time)
Memory size used (GB)Average

:
:
:
:
:
:
:

454.266
1815428.791
18425.031
1440471.743
32290.442
15181.807
4427.528<--- 15.2 TFlops

List 1: An example of MPIPROGINF output.

6
processors
grid points Tﬂops eﬃciency
4096 511 × 514 × 1538 × 2 15.2
46%3888 511 × 514 × 1538 × 2 13.8
44%3888 255 × 514 × 1538 × 2 12.1
39%2560 511 × 514 × 1538 × 2 10.350%2560 255 × 514 × 1538 × 2 9.1745%1200 255 × 514 × 1538 × 2 5.40
56%
TABLE II: Performance achieved by the yycore code—a ﬁnite diﬀerence geodynamo simulation program based on the Yin-Yang
grid— on the Earth simulator.
Paper
Shingu[16] Yokokawa[20]Sakagami[15]Komatitsch[8] Kageyama et al.
Flops/PN
26.6T/640
16.4T/512
14.9T/512
5T/243
15.2T/512
eﬃciency
65%
50%
45%
32%
46%
grid points (g.p.)
7.1 × 108
8.6 × 109
1.7 × 1010
5.5 × 109
8.1 × 108g.p./AP
1.4 × 105
2.1 × 106
4.2 × 106
2.8 × 106
2.1 × 105
Flops/g.p.
38K
19K
0.87K
0.91K
19K
Simulation kind
ﬂuid
ﬂuid
ﬂuid wave propagation
ﬂuid
Field
atmosphereturbulence inertial fusion
seismic wave
geodynamo
Methodspectral
spectralﬁnite volume spectral element ﬁnite diﬀerence
Parallelization MPI-microtask MPI-microtask HPF (ﬂat MPI)
ﬂat MPIﬂat MPITABLE III: Performances on the Earth Simulator reported at SC

To measure the performance dependency on the total number of processes, and the problem size (grid number), we
executed several runs that are summarized in Table II.
Generally, ﬂat MPI parallelization requires a larger problem size to achieve the same level of performance eﬃciency
compared to the hybrid parallelization (e.g., MPI for inter-node and microtasking for intra-node parallelization) on
the Earth Simulator [12]. Since one Earth Simulator node has 8 APs (see TABLE I), the ﬂat MPI method generates 8times as many MPI processes as hybrid parallelization. However, in our yycore code with ﬂat MPI, high performance
could be achieved with relatively low numbers of mesh size: 15.2 Tﬂops/512 PN with 511 × 514 × 1538 × 2 = 8 × 108
grid points; 12.1 Tﬂops/486 PN with 255 × 514 × 1538 × 2 = 4 × 108 grid points. This can be compared with otherﬂat MPI simulations on the Earth Simulator as reported in SC2002 and SC2003 (see Table III). Sakagami et al. [15]achieved 14.9 Tﬂops/512 PN with 1.7 × 1010 grid points for a ﬂuid simulation. (They used HPF to generate ﬂat MPI
processes.) Komatitsch et al. achieved 5 Tﬂops/243 PN with 5.5 × 109 grid points [8] for seismic wave propagation.
These results clealy indicate the high potential of the Yin-Yang grid as a base grid system on massively parallel
computers.

V.

SIMULATION RESULTS

For geodynamo study, it is necessary to follow the time development of the MHD system until the thermal convection
ﬂow and the dynamo-generated magnetic ﬁeld are both suﬃciently developed and saturated. (Initially, both theconvection energy and the magnetic energy are negligibly small.) For the case of grid size of 255×514×1538×2 = 4×108with 3888 processes, it took six hours of wall clock time until both the dynamo-generated magnetic ﬁeld and convection
ﬂow energy reached to a saturated, and balanced, level. In normalized time units, this is about 0.3% of the magnetic
free decay time. This suggests that we need to integrate more time until we observe the dynamical features of thegeodynamo such as the repeated dipole reversals [5, 11, 13].
As we described in Section III, the fundamental variables in our simulation are the magnetic vector potential A ={Ar , Aθ , Aφ }, the mass ﬂux density f = {fr , fθ , fφ }, pressure p, and the mass density ρ. It is convenient for datavisualization/analysis purpose to store the Cartesian components of the magnetic ﬁeld B = {Bx , By , Bz }, velocity
v = {vx , vy , vz }, vorticity ω = {ωx , ωy , ωz }, and temperature T . During one simulation run of 6 hours of wall clock
time, we saved the 3-dimensional data 127 times, and about 500 GB of data was generated in total.
It is known that thermal convection motion in a rapidly rotating spherical shell is organized as a set of columnar
convection cells. The number of convection columns increases and the ﬂow and the generated magnetic ﬁeld becomes7

FIG. 2: Thermal convection structure obtained by 3888 processor calculation on the Earth Simulator. The total grid size is
255(radial) × 514(latitudinal) × 1538(longitudinal) × 2(Yin and Yang). In these visualizations, grid points are reduced by afactor of 1/100. (a) Convection ﬂow in the equatorial plane viewed from the north; (b) Same data as (a) but viewed from45◦ N ; (c) Columnar convection cells viewed in the equatorial plane. Two colors indicate cyclonic and anti-cyclonic convection
columns; (d) Same data as (c), closer view.

turbulent as shown in Fig. 2, which was obtained in the simulation of 3888 processors with 255 × 514 × 1538 × 2 grid
points.

VI.

SUMMARY

In this paper, we reported the development of a new geodynamo simulation code and presented its performance on
the Earth Simulator. The base grid is the Yin-Yang grid, which is a spherical overset grid composed of two identicalcomponent grids. The achieved performance is 15.2 Tﬂops by 4096 processors with 511 × 514 × 1538 × 2 grid points

8in a spherical shell.

Acknowledgments

We would like to thank Prof. Tetsuya Sato, director of the Earth Simulator Center, for his support and encouragement. We also thank all members of the Earth Simulator Center, especially Hiromitsu Fuchigami, Shigemune
Kitawaki, Hitoshi Murai, Satoru Okura, and Hitoshi Uehara for their technical support and discussion. We also
thank Dr. Daniel S. Katz for polishing the language in the manuscript.

[1] G. Chesshire and W. D. Henshaw. Composite overlapping meshes for the solution of partial diﬀerential equations. J.
Comput. Phys., 90:1–64, 1990.
[2] Shinichi Habata, Mitsuo Yokokawa, and Shigemune Kitawaki. The Earth Simulator System. NEC Res. & Develop.,44(1):21–26, 2003.
[3] Keisuke Hirai, Keiko Takahashi, Hidenori Aiki, Koji Goto, and Kunihiko Watanabe. Development of non-hydrostatic ocean
circulation simulation code on the Earth Simulator. Proc. of The 2004 workshop on the solution of partial diﬀerential
equations on the sphere, page 65, 2004.
[4] A. Kageyama, T. Sato, K. Watanabe, R. Horiuchi, T. Hayashi, Y. Todo, T.H. Watanabe, and H. Takamaru. Computer
simulation of a magnetohydrodynamic dynamo. II. Phys. Plasmas, 2:1421–1431, 1995.
[5] Akira Kageyama, Marcia M. Ochi, and Tetsuya Sato. Flip-ﬂop transitions of the magnetic intensity and polarity reversals
in the magnetohydrodynamic dynamo. Phy. Rev. Lett., 82:5409–5412, 1999.
[6] Akira Kageyama and Tetsuya Sato. Generation mechanism of a dipole ﬁeld by a magnetohydrodynamic dynamo. Phys.
Rev. E, 55:4617–4626, 1997.
[7] Akira Kageyama and Tetsuya Sato. The “Yin-Yang Grid”: An overset grid in spherical geometry. Geochem. Geophys.
Geosys., in press; preprint: arXive:physics/0403123, 2004.
[8] Dimitri Komatitsch, Seiji Tsuboi, Chen Ji, and Jeroen Tromp. A 14.6 billion degrees of freedom, 5 teraﬂops, 2.5 terabyte
earthquake simulation on the Earth Simulator. Proceedings of the ACM/IEEE Supercomputing SC’2003 conference, 2003,
2003.
[9] Kenji Komine, Keiko Takahashi, and Kunihiko Watanabe. Development of a global non-hydrostatic simulation code using
yin-yang grid system. Proc. of The 2004 workshop on the solution of partial diﬀerential equations on the sphere, pages
67–69, 2004.
[10] Masaru Kono and Paul H. Roberts. Recent geodynamo simulations and observations of the geomagnetic ﬁeld. Rev.
Geophys., 40:1013, 2002.
[11] Jinghong Li, Tetsuya Sato, and Akira Kageyama. Repeated and sudden reversals of the dipole ﬁeld generated by a sphericaldynamo action. Science, 295:1887–1890, 2002.
[12] Kengo Nakajima. Three-level hybrid vs. ﬂat MPI on the Earth Simulator: Parallel iterative solvers for unstructured grids
on GeoFEM platform. RIST/TOKYO GeoFEM Report, 2002-007, 2002.
[13] Marcia M. Ochi, Akira Kageyama, and Tetsuya Sato. Dipole and octapole ﬁeld reversals in a rotating spherical shell:
Magnetohydrodynamic dynamo simulation. Physics of Plasmas, 6:777–787, 1999.
[14] Mitsuru Ohdaira, Keiko Takahashi, and Kunihiko Watanabe. Validation for the solution of shallow water equations in
spherical geometry with overset grid system. Proc. of The 2004 workshop on the solution of partial diﬀerential equations
on the sphere, page 71, 2004.
[15] Hitoshi Sakagami, Hitoshi Murai, Yoshiki Seo, and Mitsuo Yokokawa. 14.9 TFLOPS three-dimensional ﬂuid simulation
for fusion science with HPF on the Earth Simulator. Proceedings of the ACM/IEEE Supercomputing SC’2002 conference,
2002.
[16] Satoru Shingu, Hiroshi Takahara, Hiromitsu Fuchigami, Masayuki Yamada, Yoshinori Tsuda, Wataru Ohfuchi, Yuji Sasaki,
Kazuo Kobayashi, Takashi Hagiwara, Shin ichi Habata, Mitsuo Yokokawa, Hiroyuki Itoh, and Kiyoshi Otsuka. A 26.58
Tﬂops global atmospheric simulation with the spectral transform method on the Earth Simulator. Proceedings of the
ACM/IEEE Supercomputing SC’2002 conference, 2002.
[17] Joseph L. Steger, F. Carroll Dougherty, and John A. Benek. A Chimera grid scheme. Advances in Grid Generation, editedby K.N. Ghia and U. Ghia:59–69, 1983.
[18] Keiko Takahashi et al. Development of nonhydrostatic coupled ocean-atmosphere simulation code on the earth simulator.
Proc.of 7th International Conference on High Performance Computing and Grid in Asia Paciﬁc Region, Omiya, Japan,
2004.
[19] Peng Xindong, Feng Xiao, Keiko Takahashi, and Takashi Yabe. Conservative CIP transport in meteorological model.
JSME International Journal, to be published, 2004.
[20] Mitsuo Yokokawa, Ken’ichi Itakura, Atsuya Uno, Takashi Ishihara, and Yukio Kaneda. 16.4-tﬂops direct numerical simulation of turbulence by a Fourier spectral method on the earth simulator. Proceedings of the ACM/IEEE Supercomputing
SC’2002 conference, 2002.

9[21] Masaki Yoshida and Akira Kageyama. Application of the Yin-Yang grid to a thermal convection of a Boussinesq ﬂuid with
inﬁnite Prandtl number in a three-dimensional spherical shell. Geophys. Res. Lett., 31:doi:10.1029/2004GL019970, 2004.