Journal of the Earth Simulator, Volume 4, December 2005, 2–10

ACuTEMan: A multigrid-based mantle convection simulation codeand its optimization to the Earth Simulator
Masanori Kameyama
The Earth Simulator Center, Japan Agency for Marine-Earth Science and Technology,
3173-25 Showa-machi, Kanazawa, Yokohama, 236-0001, Japan
Email: kameyama@jamstec.go.jp(Received April 27, 2005; Revised manuscript accepted June 6, 2005)
Abstract We report the current status of our numerical simulation code named “ACuTEMan” for large-scale
mantle convection problems on the Earth Simulator. The ACuTEMan code comes out from a newly-developed “ACuTE” algorithm, which solves the flow field of mantle convection in combination with the multigrid method. This algorithm iteratively solves the steady-state equations for conservation of mass andmomentum for highly viscous and incompressible fluids, and is proved to be suitable for mantle convection
problems with excellent vector and parallel capability. In addition, a careful adaptation of the multigrid procedure enabled us to obtain a moderate parallel efficiency using up to 64 processor nodes of the Earth
Simulator. In this paper we present the optimization and the performance of our method applied to mantle
convection problems in a three-dimensional rectangular domain, with special emphasis on the enhancementof parallel efficiency of multigrid procedures by an “agglomeration” technique.
Keywords: Mantle Convection, Variable Viscosity, Multigrid Method, Parallelization, Agglomeration

1. Introduction
The convective motion in the Earth’s mantle is the ultimate origin of the geological and geophysical phenomena
observed at the Earth’s surface, such as seismicity, volcanism and plate tectonics [1, 2]. A major tool for understanding the mantle convection is numerical analysis.
Many numerical codes have been developed since thepioneering works in 1960’s and 1970’s [3, 4], and they
have been playing an important role in the study of mantle convection (see [5] for a review). Now that the Earth
Simulator is open to the studies on solid earth sciences,
further progress is expected in numerical models of mantle convection in order to deepen our understanding of the
dynamics of the Earth’s interior.
In reality, however, numerical researchers of mantle
convection seem to be suffering from the difficulties in
adapting their codes to the Earth Simulator. The difficulties mainly lie in the solution procedures of flow fields ofmantle convection. Because of the extremely high viscosity of mantle materials (~ 1022 Pa s) [1, 5], the flow in the
mantle is described by a steady-state Stokes flow balancing among the buoyancy force, pressure gradient and viscous resistance. In addition, their viscosity varies by several orders of magnitude depending on temperature, pressure, and stress [6, 7]. Taken together with the assumption of incompressibility, one needs to solve ill-condi-

2tioned elliptic differential equations for velocity and pressure at every timestep. In order to promote the mantle
convection studies in the era of Earth Simulator, it is very
important to develop efficient numerical techniques thatcan deal with the steady-state flow of highly viscous andincompressible fluids with a strongly variable viscosity.
To overcome these difficulties, we developed a new
simulation code named “ACuTEMan” for mantle convection problems on the Earth Simulator. A major innovationin our code is a solution algorithm for the flow fields ofmantle convection. This algorithm makes good use of themultigrid method [8, 9, 10] together with a new smoothing algorithm named “ACuTE” [11]. The ACuTE algorithm is proved to be suitable for mantle convection problems with an excellent vector and parallel capability.
Moreover, through a further optimization of multigrid
procedure, we obtained a moderate overall parallel efficiency with our code. In the following part of this paper,
we will introduce the solution techniques employed inour code. In particular, we will describe the solution algorithm of flow fields in detail, with special emphasis on
the optimization of multigrid method.

2. Overview of “ACuTEMan”
Our “ACuTEMan” code is basically designed for thermal convection of a highly viscous and incompressible

M. Kameyama

fluid with a strongly variable Newtonian viscosity in athree-dimensional Cartesian geometry. The nondimensional forms of the fundamental equations are (for example [12, 13]);(2)and conduction. The velocity field v at old time is used in
evaluating the advection term. The discretized equation for
T can be solved in a straightforward manner by, for example, a fully explicit, a fully implicit scheme, or their combination (Crank-Nicholson scheme). In the previous paper[11], we had already shown the results calculated by bothof a fully explicit scheme and a fully implicit scheme, andverified the reliability of our numerical results.

(3)

2.2 Multigrid Solver for Flow Field(1)

This part solves the equations (2) and (3) simultaneously for the velocity v and pressure p fields at new timeusing the new temperature T updated by (1). Here weassume, for simplicity, the viscosity
at new time isalready known by some means and kept unchanged whilethe velocity and pressure are computed. These assumptions yield a set of linear elliptic equations for velocity
and pressure, which is solved by a multigrid method.
The multigrid is known as the optimal method for anumerical solution of large-scale elliptic equations. The
basic idea of the multigrid is to solve the elliptic problem
using a hierarchy of grids with different numbers of grid
points (see Figure 1). Suppose an iterative procedure forsolving the discretized problem on the finest grid.
Assuming that the approximate solutions are already calculated on coarser grids, the iterative procedure on thefinest grid can be started from an initial guess provided byan interpolation of the approximates on the coarser grids.
Since the initial guess defined thus is most likely to be agood approximation, the computational costs of the iterative procedure is expected to be significantly saved.
However, in most of practical multigrid methods, thecomputations on different grids are intended to reduce (orsmooth) the components of errors in the approximate solutions with different spatial wavelengths. The smoothing
calculations on finer grids are responsible for reducing thecomponents with shorter wavelengths, while those on

where superscript t means transpose, v is fluid velocity, ppressure,viscosity, T temperature, Ra the Rayleighnumber, the unit vector in z-direction, and q is the internal heating rate. We assumed that z-axis is the verticalaxis pointing upward. Boussinesq approximation is
employed in the energy equation (1) and, hence, theeffects of adiabatic and viscous heating are ignored. The
temporal evolution of thermal convection is calculated
separately in two parts. In the first part, we solve theenergy equation (1), and update the temperature field. In
the second part, we solve the equation of motion (2)together with the equation of continuity (3), and update
the velocity v and pressure p fields.
In the following subsections, we will describe the solution methods employed in our code in turn. In particular,
we will discuss in detail the algorithm and its implementation of solution method for flow field (v and p), where
we made a distinct innovation for efficient computations
on the Earth Simulator.

2.1 Solver for Temperature Field
This part updates the temperature field T through a temporal integration of the energy equation (1). Equation (1) is
discretized by a first-order Euler method in time. An
upwind scheme, called power-law scheme [14], is used toevaluate the contributions of heat transport by advection

coarse↔

fineFig. 1 A schematic diagram of the principle of the multigrid method. By providing an initial guessfrom the interpolation of approximates on coarser grids, the computational cost for the iterative
procedure for solving the discretized problem on a particular grid can be significantly saved.

J. Earth Sim., Vol. 4, Dec. 2005, 2 – 10

3

Multigrid-based mantle convection simulation code “ACuTEMan”

coarser grids are for the components with longer wavelengths. By appropriately combining the approximates on
different grids, the desired solution is obtained on thefinest grid. A more detailed description of the multigrid
methods can be found in, for example, [8, 9, 10].
The multigrid procedures are, in general, composed ofmultigrid components (including smoothers, inter-gridtransfers, and derivation of coarse-grid equations) andtheir hierarchical combination into multigrid cycles overthe entire grid levels. In the rest of this subsection, we will
introduce them in turn. In particular, we will discuss indetail the smoothing operator and the multigrid cycles,
where we made special adaptations to the Earth Simulator.
2.2.1 Smoothing algorithm “ACuTE”
In this subsection, we briey introduce a newly-developed smoothing algorithm named “ACuTE”. The namecomes from an acronym of “Artificial Compressibility
using local Time Evolution”. The detail of the algorithm
can be found in [11].
The ACuTE algorithm comes from an extension of theartificial compressibility method [15], which is often used
to solve steady-state flow of incompressible fluid withhigh Reynolds number, to the problems of highly viscous
fluid with variable viscosity. In this algorithm, instead of
directly handling (2) and (3), we consider the following
set of “pseudo-evolutionary” equations for v and p:(4)(5)
Here is analogous to time, and M and K are positive
parameters analogous to density and compressibility,
respectively. The pseudo-time integration is repeateduntil a steady solution is achieved. It is obvious that thesteady solution of (4) and (5) satisfies (2) and (3), andthat it does not depend on the choice of M or K. In addition, we employ the local time-stepping technique inorder to accelerate the convergence to a steady state in
the presence of spatial variations in . In the present
algorithm, this technique is taken into account by varying
M and K depending on . As had been demonstrated in
the earlier paper [11], the convergence becomes sufficiently fast even for the cases with variable viscosity, by
assuming M ∝ and K ∝ −1. Indeed, as will be shown
later (see (6)), the convergence behavior of the ACuTE
algorithm is mainly controlled by the distribution of viscosity . This also implies that the convergence behavior
is not, in principle, affected by other parameters such as
the Rayleigh number Ra or the internal heating rate q.
We note that the ACuTE algorithm is suitable for large-

4

scale numerical simulations on massively vector-parallel
supercomputers such as the Earth Simulator. The pseudotemporal evolution of (4) and (5) can be naturally vectorized and parallelized. In particular, its vectorization andparallelization are quite straightforward if these equations
are discretized by a fully explicit scheme in the direction ofpseudo-time . In addition, this algorithm can be regarded
as a kind of iterative solvers of linear simultaneous equations. This implies that this algorithm can be directly used
as a smoother of multigrid method in place of basic iterative methods (such as Jacobi or Gauss-Seidel methods).
We also note that we must carefully choose the costs ofsmoothing calculations by (4) and (5), in particular when
the spatial variation in viscosity is involved. To see thismore clearly, we consider the Courant-Friedrichs-Lewy(CFL) criterion for (4) and (5). From these equations weobtain the pseudo-evolutionary equation of · v as,(6)(In (6) we assumed M ∝ and K ∝ −1.) The first termin the right-hand side of (6) represents the effect of“pseudo-sound wave” which propagates at a uniform rate
, and the second term represents the effect ofviscous damping at a uniform diffusion coefficient. The effect of spatial variation in is included in
the third and fourth terms, and it acts as an additional“advection velocity” of. Therefore, the CFL criterion demands smaller increment of pseudo-timein proportion to the viscosity variation. In other words, thesmoothing effect of single pseudo-temporal integration
becomes less significant as the spatial variation in isstronger. A simple remedy is to increase the number ofpseudo-temporal integrations according to the magnitudeof the spatial variation in . In the present version of
ACuTEMan, the number of pseudo-time integrations Ns
for each smoothing step on the grid level is taken to be,
(7)
where
is the magnitude of viscosity variation in theentire domain. As had been demonstrated in the earlier
paper [11], we obtained a significant improvement in therobustness of our method against the viscosity variation .
2.2.2 Multigrid method and its optimization for vector-parallel environments
One of the difficulties in optimizing the multigrid pro-J. Earth Sim., Vol. 4, Dec. 2005, 2–10

M. Kameyama

grams comes from the fact that the computations on coarser grids become less efficient. As the grid becomes coarser, the size of the discretized problems becomes smaller.
The reduction of problem sizes on coarse grids, which is
the heart of multigrid methods, degenerates the computational efficiency from two aspects. First is that theamounts of arithmetic computations become smaller as
the grid becomes coarser. This may severely spoil thecapability of vector operations, since the reduction ofproblem sizes results in shorter vector lengths on coarser
grids. Second is that the relation between computation andcommunication becomes worse than on fine grids. As the
problem size becomes smaller, the relative communication
overhead (as compared to the time for arithmetic computations) increases, and may finally dominate the arithmetic
work on the very coarse grids. This can result in a significant loss of efficiency for the overall parallel computation.
Another difficulty may arise when performing parallel
computations with a large number of processor elements(PEs). If we proceed coarser and coarser grids during a
multigrid cycle, smaller and smaller numbers of meshes
are assigned to each PE. Then at a certain coarse gridlevel, the number of meshes becomes smaller than that of
PEs involved in parallel computations. In such a situation, the domain decomposition strategy becomes incompatible with the grid coarsening for very coarse grid levels. Even if the calculations are parallelized by somemeans, the relative communication overhead on thesegrids cannot be affordable.
The above consideration suggests that the key to optimizing the parallelized multigrid is to enhance the computational efficiency particularly on coarse grid levels. In
the following, we will introduce our strategies to reduce
the time spent in the execution on coarse grids.
(1) Tips for construction of Multigrid cycles
We construct the multigrid cycle as summarized in
Figure 2. In our method, the V-cycle is used with slight
modification. The benefit of V-cycle is that it visits coarser grids at fewer times than other multigrid cycles (suchas W-, F-cycles) do. In addition, we employed an “adaptive” strategy in determining which of adjacent grid levels is visited. We measure how much the residuals are
reduced after every pre-smoothing step during the Vcycle (green circles in Figure 2). If the residuals are sufficiently reduced at a certain grid level, we go up to thefiner grid level (thick green arrows in Figure 2), insteadof going down to the coarser grid level. By this configuration of multigrid cycle, the calculations on coarse grids
can be saved as much as possible.
(2) Tips for domain decomposition
The basic strategy for parallelization is, as well as inmany applications, based on the domain decomposition

EXIT

fine

grid level

pre/post-smoothing with
fixed number of ACuTE iterations

restictionof residuals

prolongationof errors

convergence checkafter smoothingsmoothing repeateduntil convergencecoarsemultigrid steps

Fig. 2 A schematic diagram of the multigrid cycle employed in this study. Shown in red are the course of calculations in a standard V-cycle of multigrid iterations. The red dots indicate the pre- and post-smoothing steps
where a fixed number of ACuTE iterations are carried out. The red downward arrows indicate the restrictionof residuals, followed by the visit to a coarser grid level. The red upward arrows indicate the prolongation oferrors, followed by the visit to a finer grid level. The blue circle indicates that the smoothing steps are repeatedat the coarsest grid level until a convergence is achieved. The green circles indicate that a convergence checkis performed after the smoothing steps. The thick green arrows indicate the courses of calculations taken when
the residual becomes sufficiently small. In these cases we skip the visits to coarser grids, and either move backto a finer grid level or stop the multigrid iterations.

J. Earth Sim., Vol. 4, Dec. 2005, 2 – 10

5

Multigrid-based mantle convection simulation code “ACuTEMan”

method with partial overlaps. By applying this strategy toall of the grid levels involved, the multigrid components(such as smoothers) can be parallelized in a straightforward manner. In the present method, the widths of theoverlap regions are taken to be 1 for all grid levels, which
are sufficient for the spatial discretization with secondorder accuracy. Here we also note that the geometrical
domains at all grid levels should be divided in the samemanner. In other words, each PE is responsible for thecomputations on the same subdomain on a sequence ofgrid levels. By assigning the same geometric positions ofdecomposed subdomains for all grid levels, we can save
communications during inter-grid transfers.
In addition, we adjusted the manner of domain decomposition on very coarse grids where the number of computational meshes becomes small compared to that of PEs
involved in parallel computations. Figure 3 summarizes
our strategy. Suppose a model problem which is solved
using Ngrid grid levels and NPE processor elements (PEs).
Let us first divide the Ngrid grid levels into two groups,
one with grid level
and the other with(see blue dashed line in Figure 3), depending on whether
or not the number of mesheson grid level issufficiently larger than NPE (NPE
512, Ngrid
8 and
4 in the present case). Next, we apply different
patterns of domain decomposition to the different groupsof grid levels. On fine grid levels(namely), the entire computational domain isdecomposed into NPE subdomains. According to the standard parallelization strategy, the computations on thesegrid levels are carried out in parallel, by assigning each
subdomain to different PEs. On coarse grid levels where
the above condition is not met (), in contrast, wegive up decomposing the entire meshes into subdomains.
Instead of using all of the available PEs in parallel, thecomputations in the entire domain are carried out by one
PE (in this case the PE with rank 0 is chosen), leaving
other PEs idle. This strategy is sometimes called“agglomeration” technique [10]. The idea of thisapproach is to “agglomerate” the computational meshes
to new “process units” and to redistribute these new units
to a subset of the active PEs.
A major benefit of the agglomeration is that this strategy improves the overall patterns of communications. As
had been already pointed out, the main cause of communication overhead was the frequent communications during the smoothing calculations on coarse grid levels. In
the “agglomerated” case, in contrast, significant portionof communication overhead can be removed, since nocommunication is required on grid levels where only one
PE is involved in computations. Of course, the agglomeration itself requires a large amount of communication
associated with the inter-grid transfers between
and(thick green arrows in Figure 3). The data

example for 1024x1024x256 mesh problem with 512 PEs(with 2x32x8 domain decomposition)mesh divisions
in each active PE

level 8(1024x1024x256)

512x32x32

level 7(512x512x128)

256x16x16

activate
512PEslevel 6
(256x256x64)
level 5(128x128x32)
level 4(64x64x16)
level 3(32x32x8)level 2
(16x16x4)
level 1(8x8x2)

128x8x8
64x4x4

restrictionof residuals

gatherto rank 0

prolongationof errors

scatterfrom rank 0

32x2x232x32x8

activate
1PE only

16x16x4
8x8x2

V-cycle step
Fig. 3 A schematic diagram of the course of multigrid iterations together with the pattern of domain
decompositions employed in this study. The red dots and red arrows indicate the course of calculations in a standard V-cycle of multigrid method. The dotted blue line indicates a threshold
grid level which divides the entire grid levels into two groups. The grid levels finer than thethreshold contain a sufficient number of meshes compared to that of available processor elements (PEs) for parallel computations, while the coarser levels do not. The calculations on thefiner grid levels are performed in parallel using all available PEs, while the calculations on thecoarser levels are done using only one PE with rank 0. The thick green arrows indicate thecommunications between all PEs and the PE with rank 0, which take place during the transitions of grids across the threshold level.

6

J. Earth Sim., Vol. 4, Dec. 2005, 2–10

M. Kameyama

must be gathered from all PEs to one PE when going
down from grid level
and, similarly, must be
scattered from one PE to all PEs when going up to. However, the increase in communication costsby agglomeration is most likely to be compensated by thedecrease coming from the removal of frequent communications during the smoothing calculations for.
Another benefit of the agglomeration is that it can
improve the computational efficiency during the calculations on. This is somewhat counterintuitive,
since the reduction in the number of active PEs nominally
increases the execution time for the same amount arithmetic works. However, this strategy can work particularly
on vector architectures. In the agglomerated case a larger
number of meshes are mapped to active PEs than in nonagglomerated case, which results in better vectorization
because of larger vector lengths. The increase in computational time by agglomeration is most likely to be compensated as a whole by the decrease coming from theenhanced vector capability for.
2.2.3 Notes for Other Multigrid Components
Other multigrid components are constructed according
to a “standard” strategy of multigrid. A linear interpolation is used in both fine-to-coarse (restriction) and coarseto-fine (prolongation) operations. The discretized equations on coarser grids are derived by directly discretizing
the differential equations (2) and (3) on the particularcoarser grids. The values of viscosity on coarser grids,
which are necessary to derive the coarse-grid equations,are calculated by an appropriate interpolation from thoseon the finer grids by one grid level.
Here we note that the “standard” treatment describedabove may be inappropriate when viscosity varies strongly in space. More precise treatments, such as operatordependent transfer and Galerkin approximation [9, 10],are expected to improve the robustness against the viscosity variation [16].

2.3 Performance of ACuTEMan: Impact ofAdjustment in Multigrid
We measured a performance of our simulation code on
the Earth Simulator. A model problem is a thermal convection of an isoviscous fluid with Rayleigh numberRa
107 in a three-dimensional rectangular domain of
aspect ratio 4 × 4 × 1. The computational domain isdivided into equally-spaced 1024 × 1024 × 256 meshdivisions, and decomposed into 2 × 32 × 8 subdomains in
x-, y- and z-directions, respectively. The calculation isparallelized by a so-called “flat MPI” approach, and is
carried out using 64 processor nodes (i.e. 512 PEs) of the
Earth Simulator. In order to calculate the temporal evolution of thermal convection, we solve the energy equation(1) together with (2) and (3). We performed one thousand
temporal integrations, starting from a purely conductive
state plus small sinusoidal temperature perturbation. At
each timestep, the equations for v and p are solved by themultigrid method summarized in Figure 3, using
grid levels and. The mesh spacing is
successively doubled as the grid level decreases, andFig. 4 A summary of the overall performance of the ACuTEMan code. The values are reported by
MPIPROGINF available on the Earth Simulator.

J. Earth Sim., Vol. 4, Dec. 2005, 2 – 10

7

Multigrid-based mantle convection simulation code “ACuTEMan”

the mesh spacing is 1/2 on the coarsest grid level.
Figure 4 summarizes the overall performance of our code.
In this case we achieved 765 Gflops using 512 PEs, and
18% of peak performance. We also found that our code
has a moderate vectorization capability, with average
vector length 194.6 and vector operation ratio 98.6%.
To demonstrate the impact of our multigrid implementation in detail, we show in Figure 5 the plots of the execution time
per one smoothing step against gridlevel . The values are separately plotted forand(the value foris not shown, because the gridlevel
was not visited at all in this calculation). Also
shown by a thick black line in the figure is a theoretical
relation ofwith , estimated solely from thechanges in problem sizes. By comparing
the actual relations between
and with the theoretical one, we can estimate the influence of communication
overhead. Indeed, the slope of the plot ofagainstis close to that of theoretical relation forwhere nocommunication takes place during the smoothing steps.
For, in contrast, the slope significantly deviatesfrom the theoretical one as decreases. This implies thatthe influence of communication overhead becomes severer as the amount of arithmetic work decreases. In particu--1

execution time per smoothing step (sec)

10

-2

10

theoretical
relationslar, the influence is quite severe for.
From the comparison of
betweenand, we can see the impact of the agglomeration technique on the efficiency of smoothing calculations on coarse grid levels. Let us first estimate a theoretical difference inbetween4 and 3. The
amount of arithmetic works for
3 is 8 times smaller
than that for
4 owing to the reduction of problem
sizes. On the other hand, the number of active PEs for
3 is 512 times smaller than that for
4. Taken
together,for
3 is expected to be 64 times
larger than that for4. Surprisingly, however, Figure
5 shows that
increases only slightly by changing
from4 to 3. This discrepancy mainly comes from thedifference in the influences of communication overheads
between4 and 3. The smoothing steps for
3 arefree from communications, while those for
4 are significantly influenced by the communication overheads
due to the frequent communications among the active NPE(=512) PEs. The decrease in execution time by theremoval of communication dominates the potential
increase by idling the majority of PEs.
In order to see the overall parallelization efficiency in
our multigrid method, we show in Figure 6 the ratio of
total execution time of smoothing steps for all grid levels
to that of the entire smoothing calculations over one thousand temporal integrations. The portion of execution
times spent in the smoothing steps foris indicated by shading. Figure 6 clearly shows that the overall parallelization efficiency is not severely spoiled by idling thegrid level 2
68.18 sec

-3

10

-4

10

grid level 3
381.2 sec

single execution level

-5

10

7

6

5
4
grid level3

2

1
coarseFig. 5 Plots of execution times of one smoothing calculation
against the grid level obtained for a calculation of 1024 × 1024 × 256 mesh system using 512
processor elements (PEs) together with 8 grid levels and. Shown in red are the values of the finer gridlevelswhere the smoothing calculations arecarried out in parallel using all available PEs, while
shown in blue are the values ofwhere the calculations are done using one PE. We also show by athick black line, for comparison, a theoretical relations
between the execution time and grid level estimated simply from the changes in problem sizes.

grid level 4
172.8 sec

execution time forsmoothing steps
over 1000 timesteps
1923.8 sec

gr
10 id l6. eve0gr
se l 5
i
d
c
10 l
e
8. ve
9
se l 6
c

8
fine

8

grid level 1
0.000 sec

parallel execution level

grid level 8835.2 sec

grid level 7
251.6 sec

Fig. 6 A graph showing the execution times spent in thesmoothing steps for different grid levels and their ratios
to the total execution times. The shaded portion indicatesthe execution times where the smoothing calculations arecarried out using one processor element (PE) leaving
other PEs idle.

J. Earth Sim., Vol. 4, Dec. 2005, 2–10

M. Kameyama

majority of available PEs for. The smoothing
steps for these grid levels take only about one-fourth of
the total execution time, even though only one PE is
responsible for these calculations. This is a direct consequence that these calculations are executed very efficiently owing to the removal of communication (see also
Figure 5). We thus conclude that the agglomeration technique is a very effective strategy in enhancing parallelization efficiency of the multigrid method.

3. Discussion and concluding remarks
In this paper we introduced our numerical simulation
code named “ACuTEMan” for solving mantle convection
problems with strongly variable viscosity in a threedimensional rectangular box. The essence of our code liesin the improvement of the solution method for the flow
field of mantle convection, which is typically the mosttime-consuming part. In this part we employed a multigrid method, which is theoretically optimal for a solution
method for large-scale elliptic equations, in combination
with a special optimization for the Earth Simulator. We
developed a smoothing algorithm named “ACuTE” forthis problem [11], and demonstrated that this algorithm issuitable for vector and parallel computations. By optimizing the multigrid procedures using the agglomeration
technique, we successfully obtained sufficient vector andparallel efficiency of our code using up to 64 processor
nodes (i.e. 512 PEs) of the Earth Simulator. We are still
continuing to develop our code toward a more “realistic
simulation” of mantle dynamics, by incorporating theeffects of solid-state phase transitions and the spatial variations in physical properties of mantle materials [17] (seeFigure 7), and by applying our algorithm to a convection
in a three-dimensional spherical geometry [18, 19].
The performance measurement described in Section 2.3
suggest that the agglomeration technique is quite efficientin optimizing the multigrid program. In this study the
agglomeration technique is taken into account by activating only one processor element (PE) during the calculations on the grid level
while leaving other PEs
idle. One must note here that an appropriate choice ofis crucial for the overall performance of parallelized
multigrid. Suppose a case, for example, where the calculation presented in Figure 6 were performed with.
Since the problem size for
4 is eight times larger thanthat for3, it would take about (381.2 × 8 ) 3000
seconds in the calculations for4, which would haveresulted in a significant increase in execution time.
Unfortunately, however, it is impossible to derive a “golden rule” for the choice of. Since the optimalischosen to minimize the sum of the execution times forarithmetic work and those for communication, it depends,
in principle, on the specific computers as well as on thesmoothing algorithms. It is therefore necessary to findby a trial-and-error manner for specific problems.
We believe that the solution technique for flow fields
described here has a great impact on many research fields
other than mantle convection studies. In the fields ofearth sciences, efficient solvers of large-scale elliptic
equations are required in many kinds of applications,such as incompressible fluid dynamics, potential problems, and (quasi-)static viscoelastic analysis. We dedicate
this paper to Earth Simulator users who are seriously
annoyed by the solution procedure of large-scale elliptic

Fig. 7 Example of numerical simulations performed by ACuTEMan. A thermal convection in a three-dimensional rectangular box of 3000 km height and aspect ratio of 6 × 6 × 1 are calculated with 512 × 512 × 128
mesh divisions. Temperature- and depth-dependence of viscosity and temperature-dependence of thermal
diffusivity are included. An endothermic (having negative Clapeyron slope) and exothermic (positive
slope) phase transitions are imposed at the depth of around 660 km and 2800 km, respectively. Shown in
the figure are the isosurfaces of the deviation of nondimensional temperature from its horizontal average
T = −0.05 in a subdomain of 3 × 3 × 1. One temporal integration step takes about 3.4 seconds for the
calculation with 16 processor nodes (i.e. 128 PEs) of the Earth Simulator. See [17] for more details.

J. Earth Sim., Vol. 4, Dec. 2005, 2 – 10

9Multigrid-based mantle convection simulation code “ACuTEMan”

equations. Up to now, an excellent computational performance on the Earth Simulator has been reported only
on the fields where a solution of elliptic equations is not
involved [20, 21, 22, 23, 24]. We hope that our experience presented in this paper can help to broaden the horizon of a highly efficient numerical simulation of elliptic
problems on the Earth Simulator.

Acknowledgement
I thank Akira Kageyama and Tetsuya Sato for fruitful
discussion and encouragement, and Yoshinori Tsuda andKen’ichi Itakura for technical assistance. I also thankProf. Takesi Yukutake for careful and thoughtful review.
All of the calculations presented in this paper were doneby the Earth Simulator at Japan Agency for Marine-Earth
Science and Technology, as a part of the “Dynamics ofCore-Mantle Coupled System” project.
(This article is reviewed by Dr. Takesi Yukutake.)

References[1] D. L. Turcotte and G. Schubert, Geodynamics:
Applications of continuum physics to geological problems,
450 pp., John Wiley, New York, 1982.
[2] F. D. Stacey, Physics of the Earth, Third edition, 513 pp.,
Brookfield Press, Brisbane, Australia, 1992.
[3] D. L. Turcotte and E. R. Oxburgh, Finite amplitude convection cells and continental drift, J. Fluid Mech., vol. 28,
pp. 29 – 42, 1967.
[4] D. P. McKenzie, J. Roberts, and N. O. Weiss, Convectionin the Earth’s mantle, Tectonophysics, vol. 19, pp.
89 – 103, 1973.
[5] G. Schubert, D. L. Turcotte, and P. Olson, Mantle convection in the Earth and planets, 940 pp., Cambridge
University Press, Cambridge, 2001.
[6] J. Weertman, The creep strength of the Earth’s mantle,
Rev. Geophys., vol. 8, pp. 145 – 168, 1970.
[7] S.-I. Karato and P. Wu, Rheology of the upper mantle: A
synthesis, Science, vol. 260, pp. 771 – 778, 1993.
[8] A. Brandt, Multi-level adaptive solutions to boundaryvalue problems, Math. Computation, vol. 31,
pp. 333 – 390, 1977.
[9] P. Wesseling, An introduction to multigrid methods,
284 pp., Wiley, New York, 1992.
[10] U. Trottenberg, C. Oosterlee, and A. Schüller, Multigrid,
631 pp., Academic Press, San Diego, 2001.
[11] M. Kameyama, A. Kageyama, and T. Sato, Multigrid iterative algorithm using pseudo-compressibility for threedimensional mantle convection with strongly variable viscosity, J. Comput. Phys., vol. 206, pp. 162 – 181, 2005.
[12] H. Schmeling and W. R. Jacoby, On modelling the lithosphere in mantle convection with nonlinear rheology, J.

10Geophys., vol. 50, pp. 89 – 100, 1981.
[13] M. Ogawa, G. Schubert, and A. Zebib, Numerical simulations of three-dimensional thermal convection in a fluidwith strongly temperature-dependent viscosity, J. Fluid
Mech., vol. 233, pp. 299 – 328, 1991.
[14] S. V. Patankar, Numerical Heat Transfer and Fluid Flow,
197 pp., Hemisphere,Washington DC, 1980.
[15] A. J. Chorin, A numerical method for solving incompressible viscous flow problems, J. Comp. Phys., vol. 2, pp.
12 – 26, 1967.
[16] W.-S. Yang and J. R. Baumgardner. A matrix-dependenttransfer multigrid method for strongly variable viscosity
infinite Prandtl number thermal convection. Geophys.
Astrophys. Fluid Dyn., vol. 92, pp. 151 – 195, 2000.
[17] M. Kameyama, D. A. Yuen, and A. Kageyama, High-resolution 3-D numerical studies on the interplay between variable thermal conductivity and post-perovskite phase transition, Eos Trans. AGU, vol. 85, no. 47, Fall Meet. Suppl.,
Abstract MR23A – 0194, 2004.
[18] A. Kageyama and T. Sato, “Yin-Yang grid”: an overset gridin spherical geometry, Geochem. Geophys. Geosyst., vol. 5,
Q09005, doi:10.1029/2004GC000734, 2004.
[19] M. Yoshida and A. Kageyama, Application of the YinYang grid to a thermal convection of a Boussinesq fluidwith infinite Prandtl number in a three-dimensional spherical shell, Geophys. Res. Lett., vol. 31, L12609,
doi:10.1029/2004GL019970, 2004.
[20] S. Shingu, H. Takahara, H. Fuchigami, M. Yamada,
Y. Tsuda, W. Ohfuchi, Y. Sasaki, K. Kobayashi,
T. Hagiwara, S.-I. Habata, M. Yokokawa, H. Itoh, and
K. Ohtsuka, A 26.58 Tflops global atmospheric simulation
with the spectral transform method on the Earth Simulator,
In Proceedings of the ACM/IEEE Supercomputing SC2002
conference, 2002.
[21] H. Sakagami, H. Murai, Y. Seo, and M. Yokokawa, 14.9
TFLOPS three-dimensional fluid simulation for fusion science with HPF on the Earth Simulator, In Proceedings of
the ACM/IEEE Supercomputing SC2002 conference, 2002.
[22] M. Yokokawa, K.-I. Itakura, A. Uno, T. Ishihara, and Y.
Kaneda, 16.4-Tflops direct numerical simulation of turbulence by a Fourier spectral method on the Earth Simulator,
In Proceedings of the ACM/IEEE Supercomputing SC2002
conference, 2002.
[23] D. Komatitsch, S. Tsuboi, C. Ji, and J. Tromp, A 14.6 billion degree of freedom, 5 teraflops, 2.5 terabyte earthquake simulation on the Earth Simulator, In Proceedingsof the ACM/IEEE Supercomputing SC2003 conference,
2003.
[24] A. Kageyama, M. Kameyama, S. Fujihara, M. Yoshida,
M. Hyodo, and Y. Tsuda, A 15.2 TFlops simulation ofgeodynamo on the Earth Simulator, In Proceedings of the
ACM/IEEE Supercomputing SC2004 conference, 2004.

J. Earth Sim., Vol. 4, Dec. 2005, 2–10