数理解析研究所講究録 1441 巻 2005 年 129-140

129
ROBUST AND EFFICIENT PARALLEL ITERATIVE SOLVERS FOR ULTRALARGE SCALE STRUCTURAL ANALYSIS
吉村 忍 (Shinobu YOSHIMURA)
Institute of Environmental Studies, University of Tokyo
九州大学・知能機械システム部門 塩谷 隆二 (Ryuji SHIOYA),
荻野 正雄 (Masao OGINO)

東京大学・環境学専攻

Department of Intelligent Machinery & Systems Engineering,
Kyushu University

1. INTRODUCTION
Various general-purpose computational mechanics systems have been developed in
the last three decades to quantitatively evaluate mechanical / physical phenomena such as
deformation of solid, heat transfer, fluid flow and electromagnetics. Nowadays such systems
are regarded as infrastructural tools for the present industrialized society. The existing
systems, however, can not be used with massively parallel processors (MPPs) with the orderof 100-10000 processing elements (PEs), as they were developed for single-processor
computers. Neither can the current systems be used in heterogeneous parallel anddistributed computing environments such as the Grid. Owing to the fact, they can deal with
only mediuinscale problems with millions degrees of freedoin(DOFs) at most.
project $[1, 2]$ was one of the research projects in the
The“Computational Science & Engineering” field selected for the “Research for the Future(RFTF)” Prograinsponsored by the Japan Society for the Promotion of Science (JSPS) [3,41during 1997-2002. The project is continuously going on as an open source softwaredevelopment project. In the project we have been developing an advanced general-purpose
computational mechanics systeinnamed ADVENTURE since August 1997. The systeinwas
designed to be able to analyze a three-dimensional finite element model of arbitrary shapewith 10-100 million DOF mesh, and additionally to enable parametric and non-parametric
The first version of the ADVENTURE systeinhas been
shape optimization $[5, 6]$ .
released ffointhe project website [1] as open source software since March 2002. About
1800 registered users in academia and industries are now using the programs, while one
private company has developed and released its commercial version named
ADVENTUREcluster $[5, 6]$ .
Domain-decomposition-based parallel algorithms are implemented in pre-processes
solutions) and(domain decomposition), main processes (systeinmatrix assembling anddecomposition
post-process (visualization), respectively. Especially the hierarchical domain
the mainmethod with a preconditioned iterative solver (HDDM) [9-12] is adopted in two of
and
Solid
modules for solid analysis and thermal conduction analysis, named ADVENTURE
The employed preconditioner is the Balancing Domain
ADVENTURE Thermal.
Decomposition (BDD) type method [13-18], To efficiently solve a coarse space problemofderived ffoinequilibriuinconditions for singular problems associated with a number
parallelbased
subdornains appeared in the BDD fomulation, an incomplete factorization$\mathrm{A}\mathrm{D}\mathrm{V}\mathrm{E}\ovalbox{\tt\small REJECT}$

130
direct solver is employed. The ADVENTURE Solid has been successfully implemented on,a single PC PC clusters and massively parallel processors such as Hitachi
11, 12, 19] Recently, this solid analysis module has been implemented with minor
modification on the Earth Simulator consisting of 256 nodes, i.e. 2,048 vector-type PEs oftheoretical peak performance of 16 TFLOPS, and succeeds in solving an elastostatic problem
of a nuclear pressure vessel model of 100 million DOFs in 8.5 minutes with 5.1 TFLOPS,
which is 31.S% of the peak performance and over 80% parallel efficiency.
$\mathrm{S}\mathrm{R}8000/\mathrm{M}\mathrm{P}\mathrm{P}[2$EA4OF 4t-IIVENT
The ADVENTURE system consists of pre-, main- and post-processing modules anddesign modules that can be used in various kinds of parallel and distributed environments [1,2], The system employs a hierarchical domain decomposition method (HDDM) [9-12]
based massively parallel algorithm as one of the major solution algorithms in order to handle
a huge-scale finite element model over 10-100 million DOFs efficiently. The system
The pre-process modules
employs module-based architecture and consists of 19 modules.
include the surface patch generator which converts geometry model data into a collection of, a tetrahedral mesh generator
triangular surface patch data, named ADVENTURE$[20, 21]$ , i.e. ADVENTURE, an attachment tool of boundary conditions andmaterial properties onto the mesh, i.e. ADVENTURE BCtool, and a domain decomposer of afinite element model, i.e. ADVENTURE Metis. The kernels of the ADVENTURE Metis
are a graph partitioning tool METIS and its parallel version ParMETIS developed in the
University of Minnesota $[22, 23]$ . The main process modules, . solvers include an implicit
elastic-plastic analysis module named ADVENTURE Solid [11-12, 19] which enables
large-deformation and implicit dynamic analyses, a thermal conductive analysis module
named
ADVENTURE Thermal,a thermal-fluid analysis module named
ADVENTURE Fluid, a magnetic analysis module named ADVENTURE Magnetic [24], an
explicit impact analysis module named ADVENTURE Impact, and a rigid plastic analysis
The post process module named
module named ADVENTURE Jorge.
ADVENTURE Visual is for parallel visualization of analysis results [25]. Figure 1 shows
the configuration of the ADVENTURE modules.

2.

$\mathrm{U}1^{-}\mathrm{t}\mathrm{E}$$\mathrm{O}1^{\gamma}\mathrm{E}\mathrm{J}^{1_{\iota}}\mathrm{t}^{\gamma}\mathrm{L}\mathrm{E}\mathrm{V}\mathrm{V}$$\mathrm{S}\mathrm{V}\mathrm{S}\mathrm{I}^{\urcorner}$$-\mathrm{T}\mathrm{r}\mathrm{i}\mathrm{P}\mathrm{a}\mathrm{t}\mathrm{c}\mathrm{h}$$-\mathrm{T}\mathrm{e}\mathrm{t}\mathrm{M}\mathrm{e}\mathrm{s}\mathrm{h}$$\mathrm{i}.\mathrm{e}$Fig. 1

Configuration ofADVENTURE modules

131]’
ADVENTURE Solid-JJVfBOne of the key technologies of the ADVENTIJRE Solid is the HDDM, whichenables parallel finite element calculations on various kinds of computing environments[9-12]. Basically in the HDDM, force equivalence and continuity conditions among
subdomains are satisfied through iterative calculations such as the Conjugate Gradient (CG)method. Therefore it is indispensable to reduce the number of iterations by adopting someappropriate preconditioning technique especially for solving large-scale ill-conditioned
problems. The Neumann-Neumann algorithm (N-N) [13] is known as efficient domain
decomposition preconditioner for unstructured subdomains. However, its convergence
deteriorates with the increasing number of subdomains due to lack of a coarse space problem
which takes care of global propagation of error. The Balancing Domain Decomposition(BDD) based N-N algorithm proposed by Mandel [14] shows that the equilibrium conditions
for the singular problems on subdomains result in simple and natural construction of a coarse
space problem and that its construction is purely algebraic. The BDD has been applied tosolve various phenomena $[15, 16]$ . There are also several researches on parallelization of the
BDD and also the FETI (Finite Element Tearing and Interconnecting) [26-30], However,most problems solved there are still medium scale ones such as sub-millions to one million
DOFs. As the DOFs of the coarse space problem is directly related to the number of
subdomains, it is indispensable to consider the parallelization of the solution process of thecoarse space problem as well when solving large-scale problems. The Salinas system [31],which employed the FETI-DP method [30], is succeeded in solving large-scale problems suchas over 100 million DOF mesh of optical shutter model [32]. It shows good performance but
does not seem to include load-balancing techniques. In the present study, an incomplete
parallel direct method and the HDDM are adopted.

3.

$\mathcal{P}\mathrm{A}\mathrm{f}\mathrm{K}\mathrm{A}\mathrm{L}\mathrm{L}\Gamma_{\lrcorner}^{4}\mathrm{L}$$\mathrm{A}\mathrm{f}\lrcorner C_{\mathrm{I}}\cdot \mathrm{O}1\mathrm{t}^{\backslash }(1^{r}.\mathrm{f}\mathrm{A}$$\mathfrak{M}$Fig.2(a) 35 million DOF mesh of

ABWR model$\mathrm{L}\mathrm{E}\mathrm{M}\mathrm{E}\grave{1}^{\mathrm{t}_{\mathrm{Q}}^{1}}\mathrm{T}\mathrm{E}\neg \mathrm{D}$$\mathrm{t}\}_{\backslash }^{\mathrm{T}}$Fig.2(b) Part decomposition

vessel mode$\mathrm{o}\mathrm{f}\mathrm{A}\mathrm{B}\mathrm{W}\mathrm{P}\backslash$132.
schematic data flow among
According to the design concept of the HDDM, most computation is assigned to theChildren, while most communication occurs in between Parents and Children. Varying thenumber of Parents and Children for different kinds of parallel computers, the present
HDDM based system can easily achieve high parallel performance. In the HDDM
architecture, thanks to the dynamic load balancing technique among Child processors, highparallel performance can be achieved even in heterogeneous computer environments.
However in this mode, an amount of data communication between Child and Parent tends to$\mathrm{P}\mathrm{E}\mathrm{s}$$\mathrm{P}\#\mathrm{r}\epsilon \mathrm{n}\mathrm{t}$$\not\in---\gg\ovalbox{\tt\small REJECT}$$\mathrm{I}$Disk$<---\gg$$i$$\ll---->\ovalbox{\tt\small REJECT}\dot{}$$\mathrm{m}$Coinmunication

$<->$

File

$\mathrm{I}/\mathrm{O}$

Fig.3 (a) Schematic data flow in h-mode

Fig.3(b) Schematic data flow in p-mode

133
be large. To reduce such data communication among Children and Parents, it is useful tois called parallel
adopt static load balance. This analysis mode shown in Figure
is named
processors mode ( -mode), while the original analysis mode as shown in Figure
hierarchical processors mode (h-mode)$3(\mathrm{b})$$3(\mathrm{a})$$\mathrm{p}$Domain Decomposition (BDD)
The BDD algorithm is based on the DDM with a preconditioned iterative solver.
After eliminating interior DOFs of local subdomain matrices, the problem to be solved isreduced onto the interface DOFs of subdomains. The reduced matrix is so-called Schurcomplement, The reduced problem is also called the interface problem, and is to be solved
by a preconditioned iterative method. There are two main methods as preconditioner, i.e.
local subdomain correction and coarse grid correction in a coarse space. Main elemental
calculations appeared in the BDD algorithm are described as follows :3.2$\mathrm{B}3\mathrm{a}\mathrm{l}\mathrm{a}\mathrm{n}\mathrm{c}\mathrm{i}\mathrm{n}\mathrm{g}$is a linear algebraic system to be solved, whereassumed to be positive definite.
a)$Ku$ $=f$b)$Su_{B}=g$

is the reduced system, whereassumed to be positive definite, and$\mathrm{i}=1,\ldots,N$c)is the Schur complement,is the local Schur complement of i-th subdomain

is the matrix of the global to local DOFs corresponding to interface mapping for i-th

$R_{i}$$R_{i}R_{i}^{T}=I$

is a weighting matrix for i-th subdomain, assumed to form decomposition of unity and$D_{i}$

satisfye)$S= \sum^{\mathrm{A}^{r}}i=1R^{r}\mathrm{i}S:R_{i}$is the global stiffiess matrix,, assumed to be positive semi-definite,

subdomain, assumed to satisfyd)$S_{i}$$K$$Z_{i}$$\Sigma_{i=1}^{N}R_{i}^{T}D_{i}R_{i}=I$isthe local coarse space of i-th subdomain, that contains all potential local

singularities,f)$R_{0}$

is the weighted restriction from the global to coarse DOFs, defined by$R_{0}^{T}=[R_{1}^{T}D_{1}^{T}Z_{1},\ldots,$ $R_{N}^{T}D_{N}Z_{\lambda’}]$

g)

$P$

$P=QS$ ,

is

the

where

positive definite.

.

-orthogonal$S$$Q=R_{0}^{T}S_{0}^{-\mathrm{J}}R_{\mathrm{I}\mathrm{J}}$, andprojection$S_{0}=R_{\mathrm{o}}SR_{0}^{T}$

onto

the

coarse space, defined by

is a coarse grid operator, assumed to be134
Various domain decomposition methods contain a process of solving a reduced
system using iterative methods such as the preconditioned CG method. At each step, theDDM or the HDDM requires to solve the following auxiliary problem :$lffz=r$(1)

where $ltf$ i a symmetric positive definite matrix called preconditioner and$\mathrm{s}$

$=F+$

where
$P^{T}r=0$$T_{\mathrm{i}}$is a residual

The BDD preconditioned operator is described $[14, 16]$ by .

vector in each iterative step.
$\mathit{1}1f^{-1}S$

$r$

(2)

$(I -p^{\backslash } \int_{\backslash }\Sigma_{\mathrm{i}1}^{N}T_{i}\mathrm{X}J -P)^{T}$

is the local subdomain correction and$J-J^{\mathrm{J}}$is the coarse grid correction.

If, which means a residual vector has no components of the coarse space, Eq. (22) canbe simplified as :(3)

$l1I^{-\mathrm{I}}S=(I-P\mathrm{X}\Sigma_{i1}^{N}T_{l})$The original BDD employs the N-N type algorithm as local subdomain correction
with a two-level weighted sum of the inverses of.

$S_{*}$

matrices [14].

To calculate the inverseof them, the Moore-Penrose pseudo-inverse or some regularization is required since$S_{\mathrm{i}}$matrices are typically singular.

However the Moore-Penrose pseudo-inverse takes highcomputational cost, while the regularization is less accurate. To overcome both issues ofcomputation cost and accuracy simultaneously, we choose the diagonal scaling preconditioner

as local subdomain correction. Since the local subdomain correction is applied
subdornain-wise, its parallel algorithm is basically compatible to the HDDM.

for

3.3

$s_{i}$

lfeliz.uliolJ of BDI)The coarse grid correction for elastic problems becomes the Galerkin projection onto
a coarse space derived from rigid body mode. For parallelizing the BDD we decompose the$\mathrm{I}^{1}\mathrm{t}\mathrm{t}\mathrm{r}‘ \mathrm{B}$BDD algorithm into the following two stages :construction of the coarse grid operator$s_{0}$’,

isAt first, the coarse grid operator
applied in subdomain-based blocks, and then its process is completely parallelized
and operation of the coarse grid correction
subdomain- wise.

That is,$\mathrm{A}\tau_{0}$

$I$

$-P$

is described by:.

$\iota 5_{\mathrm{t})}$135$S_{0}=P_{\mathrm{W}}( \sum_{k=1}^{\mathit{4}\backslash }.R_{k}^{T}S_{k}R_{\hslash})P\kappa_{0}^{T}$$= \sum_{i=1}^{N}$

where$\mathrm{i}$(4)$\mathrm{I}_{jarrow 1}^{\mathit{4}\backslash ^{7}}-\sum_{\mathrm{f}=\mathrm{I}}^{\mathit{1}\mathrm{v}}lP_{\mathrm{W}i}P_{\mathrm{t}_{k}}^{T}S_{k}R_{k}F_{\mathrm{t})j}^{T}$$P_{\mathrm{t}_{0i}}$

is the weighted restriction from the global to the local coarse DOFs of subdomain, defined by$P\iota_{0i}=[Z_{i}^{T}I\mathit{1}_{i}^{T}R_{i}]$.

Figure$4(\mathrm{a})$schematically illustrates the data flow in parallel

construction of the coarse grid operator for the HDDM.
Secondly, the coarse grid correction is applied to solve a linear system equation
whose coefficient matrix is derived from the coarse grid operator. It should be noted herethat this process can not be parallelized subdomain-wise. Thus we employ the follo wing
method. The coarse grid correction is implemented in each iteration with its own right-handside vector. Here in order to reduce computation time, an LU factorized coarse matrix in thefirst iteration can be kept, and then the forward elimination and the backward substitution of
the coarse system are applied after the second iteration.
The number of DOFs of the coarse space is directly related to the number of
subdomains, and tends to become large when solving large-scale problems. To solve thecoarse space problem, we adopt an incomplete factorization based parallel direct method.
The coarse grid operator is factorized incompletely. In general, such an incompletely
factorized operator is used together with some iterative calculation to compensate theFig.4(a) Data flow of constructing coarse grid correction

Fig.4(b) Analysis flow of the BDD with imcomplete coarse grid correction

13G
incompleteness. In this system, however, the coarse space problem is solved by the
incompletely factorized operator without iterations. This incomplete factorization process
decreases computation cost for each iteration and improves parallel efficiency, but may
increase the number of iterations. Finally computation time is expected to be reduced. In
the original BDD preconditioner of Eq. (33), the coarse grid correction is implemented after
local subdomain correction in each iteration. However, in the present BDD preconditionerof Eq. (22), the coarse grid correction is applied to the CG residual vector before local
shows the analysis flow of the present parallel BDD
subdomain correction. Figurealgorithm.
$4(\mathrm{b})$4. NUMERICAL EXPERIMENTS
4.1 Elastostatic Analysis of ABWR Vessel Model with 35 Million DOFs on
This section first describes an elastostatic stress analysis for a precise model of anadvanced boiling water reactor (ABW R) vessel with a35 million DOFs unstructured mesh as. Size of fine elements placed in nozzle corners and intem al pumpshown in Figure, while that of average elements is about lOmm. As boundaryjunctions is about
conditions, the bottom surface of its skirt portion is fixed, and a static gravitational force isapplied to the vessel in the horizontal direction, imitating a seismic loading condition. Sucha complex shaped and large-scale thin structure with less constraint often results in an
ill-conditioned system matrix. Most iterative solution methods fail in attaining convergencesuccessfully
when solving such ill-conditioned problems. The
overcomes this problem owing to the employment of the BDD-based preconditioner. Figure
5 shows convergence histories of force imbalance measured at the interface of subdomai , i.e.
residual norm, plotted against the number of iterations. The calculations are performed on
consisting of 1 ,whose theoretical peak performance is
Hitachi. There are two lines in the figure. HDDM denotes the result obtained using
the HDDM with a simple diagonal-scaling preconditioner, while BDD denotes that of the
HDDM with the BDD preconditioner. As for the HDDM case, the residual norm doesn’t$\mathrm{S}\mathrm{R}8000/\mathrm{M}\mathrm{P}\mathrm{P}$$2(\mathrm{a})$$2\mathrm{m}\mathrm{m}$$\mathrm{A}\mathrm{D}\mathrm{V}\mathrm{E}\mathrm{N}\mathrm{T}\mathrm{U}\mathrm{R}\mathrm{E}_{-}\mathrm{S}\mathrm{o}1\mathrm{i}\mathrm{d}$$\mathrm{n}\mathrm{s}$$024\mathrm{P}\mathrm{E}\mathrm{s}$$\mathrm{S}\mathrm{R}8000/\mathrm{M}\mathrm{P}\mathrm{P}$$1.8\mathrm{T}\mathrm{F}\mathrm{L}\mathrm{O}\mathrm{P}\mathrm{S}$Relate residual norm
$1.\mathrm{O}\mathrm{E}+02$$1.0\mathrm{E}+\mathrm{O}1$$\mathrm{I}.\mathrm{O}\mathrm{E}+\mathrm{O}\mathrm{O}$HDDM

1.OE-OI

I.OE-02
I.OE-03

1.OE-04
1.0E-05$\mathrm{B}\mathrm{D}\mathrm{D}$1.OE-06
1,OE-07

1

1001

2001

3001

Fig.5 History of residual norm on

4001

5001

$\mathrm{S}\mathrm{R}8000/\mathrm{M}\mathrm{P}\mathrm{P}$

6001

7001

with 1024 PEs

137
decrease even below l.Oe-l though spending 5 hours for over 7000 iterations. On the other
hand, in the BDD case, the calculation is successfully converged with about 1400 iterationsin about 20 minutes. As this model is divided into about 30000 subdomains, the number of
DOFs of the coarse space becomes about 180,000. This coarse space problem is solved inonly 30 seconds with the incomplete factorization based parallel direct method.
Figure 6 shows total computational time . the number of processors employed. It
can be evaluated ffom the figure that the scalability in total computation time is less than
99.7%, referring the value for the 128 PE case. The scalability in the total computation timeis affected by the increased number of CG iterations, which is caused due to the employmentof the incomplete factorization. It is anticipated that the present system still has somepotential to attain further speed-up by improving the incomplete factorization based parallel
direct method. Figure 7 shows calculated stress distribution and enlarged deformation of theseismic loading. By solving a whole reactor
ABW vessel subjected to thevessel with a sufficiently large scale finite element mesh, various local stress concentration
regions can be clearly indicated, and precisely evaluated.
$\mathrm{v}\mathrm{s}$$\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{s}\mathrm{i}rightarrow \mathrm{s}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{c}$$\mathrm{W}\mathrm{R}$42 Elastostatk Analysis of Pressure Vessel Model with 100 Million BOF on the ES
is implemented on the Earth Simulator [33] consisting ofTheof main memory, whose theoretical peak performance is
256 nodes, i.e. 2,048 PEs with
16 TFLOPS. The second problem is an elastostatic stress analysis of a simplified pressure
vessel model with 100 million DOFs unstructured mesh. Its mesh size is listed in Table 1.
As a boundary condition, the bottom surface of the vessel is fixed, and a static gravitational
force is applied to the vessel in the horizontal direction, being similar to the previous problem.
Although we do not show convergence histories of relative residual, (1) DDM with BDD
and N-N preconditioner (denoted as BDD) and (2) HDDM with BDD and diagonal-scaling
preconditioner (denoted as BDD-DIAG) demonstrate excellent performance in convergence.
By considering the performance results, it is concluded here that the diagonal scaling issufficient as local subdomain correction in the BDD method. The amalysis model is divided$\mathrm{A}\mathrm{D}\mathrm{V}\mathrm{E}\mathrm{N}\mathrm{T}\mathrm{U}\mathrm{R}\mathrm{E}_{-}\mathrm{S}\mathrm{o}1\mathrm{i}\mathrm{d}$$4\mathrm{T}\mathrm{B}$

3500
3000

2500
$\overline{\underline{\S \mathrm{f}\mathrm{f}\mathrm{i}}}$2000$.\succ\underline{\not\in}1500$

100

soo
0

128 256 384 512 768 896 1024$\#$

processors

Fig.6 Scalability in total computation time

Fig.7 Stress distribution and deformation$(\mathrm{x}4,000)$of ABWR vessel subjectedquasi-staic seismic loadin

138
into 34816 subdomains and then the number of DOFs of its coarse space is 208,896. The
LU factorization of the coarse grid operator is calculated in only 20 seconds. As the result,the present system successffilly achieved 5.1TFLOPS, which is 31.8% of the peak
performance. The calculation time is only 8.5 minutes. Parallel ratio over 99.9% isachieved, and then parallel efficiency exceeds 80% not only for computation time periteration but also for total computation time.
Table 1: Mesh size for a simplified vessel model
25,084,456
Number of elements
34,772,634
umber of nodes
104,195,500
Total degrees of freedom$\mathrm{N}$In the DDM algorithm, the number of subdomains employed influences performance, we have the quasi-optimu
of the calculation to some extent. In the case of
number of DOFs per subdomain from our experience, that is 600. Considering that the ES, the larger
has vector-type processors and larger memory space compared with$\mathrm{S}\mathrm{R}8000/\mathrm{M}\mathrm{P}\mathrm{P}$$\mathrm{m}$$\mathrm{S}\mathrm{R}8000/\mathrm{M}\mathrm{P}\mathrm{P}$numbers of subdomains are tested. It is estimated from some preliminary results that the.
semi-optimum number of subdomains for the ES is about 3 ,$000\mathrm{D}\mathrm{O}\mathrm{F}\mathrm{s}/\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{d}\mathrm{o}\mathrm{m}\mathrm{a}\acute{\iota}\mathrm{n}$5. CONCLUSIONS
We have been developing an advanced general-purpose finite element analysis
system, named ADENTURE, which is designed to be able to analyze a model of arbitrary
shape with a 10-100 million DOF mesh. After the overview of the ADVENTURE system isbriefly reviewed, the parallel solution algorithm, i.e. the Hierarchical Domain Decomposition
Method with Balancing Domain Decomposition based preconditioner is explained, The, PC clusters and
ADVENTURE Solid has been successfully implemented on a single
and the Earth Simulator. The
massively parallel processors such as Hitachi
ITFLOPS,
31.8%
of the peak performance fo
whichispresent system successfully achieved 5.
. The calculation time of
the Earth Simulator consisting of 256 nodes, i.e. 2,048
elastostatic analysis with 100 million DOF mesh is only 8.5 minutes. Parallel ratio over
99.9% is achieved, and then parallel efficiency exceeds 80% not only for computation time
per iteration but also for total computation time. It is clearly concluded from those resultsthat the ADVEN TURE system on the Earth Simulator can be used for virtual mockup tests oflarge-scale and complex artifacts such as nuclear pressure vessels.
$\mathrm{P}\mathrm{C}$$\mathrm{S}\mathrm{R}8000/\mathrm{M}\mathrm{P}\mathrm{P}$$\mathrm{P}\mathrm{E}\mathrm{s}$Acknowledgments
The authors would like to thank the Earth Simulator Center for giving us theopportunity to use the Earth Simulator. They would also like to thank Mr. T. Fukuda at
Tokyo Electric Power Company, Incorp. for providing us an image model of ABWR vessel.
This work is performed as a part of the ADVENTURE project [1]. The authors would also
like to thank all the members of the ADVENTURE project

139
References-tokyo.ac.jp[1][2] S. Yoshimura, R. Shioya, H. Noguchi and T. Miyamura, Advanced general-purpose$\mathrm{h}\mathrm{t}\mathrm{t}\mathrm{p};//\mathrm{a}\mathrm{d}\mathrm{v}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{e}.\mathrm{q}.\mathrm{t}\mathrm{u}$computational mechanics system for large-scale analysis and design, Journal ofComputational and Applied Mathematics, 49 (2002) 279-296.
[3] Report on Computational Science and Engieering, JSPS-RFTF Program 2001-2002,(2002).

.is. . -tokyo.ac.jp/cse/index-e.
[4][5] S. Yoshimura, N. Oshio, H. Kawai, H. Azegami and K. Takeuchi, Large-scale parallel$\mathrm{h}\mathrm{t}\mathrm{m}1/’$$\mathrm{h}\mathrm{t}\mathrm{t}\mathrm{p}:,//\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{t}\mathrm{o}\mathrm{n}$$\mathrm{s}$$\mathrm{u}$shape optimization using hierarchical domain decomposition with traction method, Journal of
the Simulation Technology, 20-4 (2001) 296-302 (in Japanese).
YoshimuraS. Dulikravich,[6] B. H. Dennis, I. N. Egorov, H. Sobieczky,
blades,inturbinepassages
Thermoelasticity optimization of 3-D serpentine coolingInternational Journal of Turbo & Jet-Engines, 21-1 (2004) 57-68.
Yagawa, Analysis of stress intensity factor of[7] M. Suzuki, H. Akiba, S. Yoshimura and,piping using large scale analysis code ADVentureCluster, Trans. 16th SMiRT,
Washington D.C., CD-ROM, (2001).
[8] M. Suzuki, T. Ohyama, H. Akiba, H. Noguchi and S. Yoshimura, Development of fast androbust parallel coarse-grid based CG solver for large scale finite element analyses, Trans.
JSME, 68A-67l (2002) 1010-1017 (in Japanese).
[9] Yagawa and Shioya, Parallel finite elements on a massively parallel computer withdomain decomposition, Computing Systems in Engineering , 4 (1994) 495-503.
[10] Shioya and Yagawa, Parallel finite elements of ten-million DOFs based on domain
decomposition method, WCCM IV Computational Mechanics -New Trends and ApplicationsIV 11 (1998) 1-12.
[11] R. Shioya and G. Yagawa, Parallel finite element analysis of 100 million DOFs based on
the hierarchical domain decomposition method, Trans. JSCES, 3 (2001) 201-206 (in$\mathrm{S}\underline{.}$$\mathrm{G}$$\mathrm{G}$$\mathrm{G}05/5$$\mathrm{K}$Japanese).
[12] T. Miyamura, H. Noguchi, R. Shioya, S. Yoshimura and G. Yagawa, Elastic-pla stic

analysis of nuclear structures with millions of DOFs using the hierarchical domain
decomposition method, Nuclear Engineering & Design, 212 (2002) 335-355.
[13] Y. H. DeRoeck and P. LeTallec, Analysis and test of a local domain decomposition
preconditioner, 4th International Symposium on Domain Decomposition Methods (1991)
112-128.
[14] J. Mandel, Balancing domain decomposition, Communications on Numerical Methods inEngineering , 9 (1993) 233-241.
[15] P. LeTallec and M. Vidrascu, Generalized Neumann-Neumann preconditioned for
iterative substructuring, 9th International Symposium on Domain Decomposition Methods,(1994) 413-425.
[16] P. LeTallec, J. Mandel and M. Vidrascu, A Neumann-Neumann domain decomposition
for solving plate and shell problems, SIAM J. Number. Math., 35 (1997)

algorithm

836-867.

Ogino,
using a balancing domain decomposition method, Trans.

[17] R. Shioya, H. Kanayama, D. Tagami and$\mathrm{M}_{\wedge}$

large scale structural analysis
JSCES, 2 (2000) 139-144 (in$3\mathrm{D}$140Japanese).
[18] R. Shioya, M. Ogino, H. Kanayama and D.Tagami, Large scale finite element analysis
with a balancing domain decomposition method, Key Engineering Materials, 243-244 (2003)

21-26.

[19] T. Miyamura and S. Yoshimura, Parallel stress analyses of ancient architecture Pantheonon PC cluster, Trans. Architectural Institute of Japan, 55 (2001) 95-102 (in Japanese).

[20] G. Yagawa, S. Yoshimura and K. Nakao, Automatic mesh generation of complex

geometries based on fuzzy knowledge processing and computational geometry, Integrated
Computer-Aided Engineering 2 (1995) 265-280.
Yagawa and H. Akiba Parallel automatic mesh generation of[21] S. Yoshimura, H. Nitta,
nuclear structures with ten-million nodes, Trans. 15th SMiRT, Seoul, II (1999) 21-28.
[22]
Karypis and V. Kumar, Multilevel k-way partitioning scheme for irregular graphs,
Technical Report TR 95-064, Department of Computer Science, University of Minnesota,$\mathrm{G}$$\mathrm{G}$(1995).
[23] G. Karypis and V. Kumar, Parallel multilevel k-way partitioning scheme for irregulargraphs, Technical Report TR 96-036, Department of Computer Science, University of

Minnesota, (1996).
eddy[24] H. Kanayama, R. Shioya, D. Tagami and M. Saito, Numerical analysis ofcurrent problems by the hierarchical domain decomposition method, Trans. JSCES, 3 (2001)
151-156 (in Japanese).
Yagawa, Parallel visualization of[25] S. Shoui, S. Yoshimura, H. Akiba, T. Ohyama andfinite element solutions with ten million DOFs using PC cluster, Proceedings of European
Congress on Computational Methods in Science and Engineering (ECCOMAS2000),
Balcelona, CD-ROM, (2000)[26] M. Vidrascu, Remarks on the implementation of the generalized neumann-neumannalgorithm, 11th International Conference on Domain Decomposition Methods, (1998)$3\mathrm{D}$$\mathrm{G}$485-493.
and P. Plechac, PARASOL. An Integrated
programming environment for parallel sparse matrix solvers, High-Performance Computing,(1999) 79-90.
[28] F.-X. Roux and C. Farhat, Parallel implementation of the two-level FETI method, 9thInternational Conference on Domain Decomposition Methods, (1997)[29] J. Mandel, R. Tezaur and C. Farhat, A scalable substructuring method by lagrange
multipliers for plate bending problems, SIAM Journal of Numerical Analysis, 36 (1999)[27] P. R. Amestoy, I. S. Duff, J. -Y$\mathrm{I}/\mathrm{E}\mathrm{x}\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\mathrm{e}\mathrm{n}\mathrm{t}$1370-1391,[30] M. Lesoinne and K. Pierson, FETI-DP :An efficient, scalable and unified dual-primal
FETI method, 12th International Conference on Domain Decomposition Methods, (1999)$421\mathrm{A}28$.

[31].sandia.
[32] M. Bhardwaj, K. Pierson, G. Reese, T. Walsh, D. Day, K. Alvin, J. Peery, C. Farhat and$\mathrm{h}\mathrm{t}\mathrm{t}\mathrm{p}://\mathrm{e}\mathrm{n}\mathrm{d}\mathrm{o}$$\mathrm{g}\mathrm{o}\mathrm{v}/9234/\mathrm{s}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{s}$M. Lesoinne, Salinas : A scalable software for high-performance struc 軸面 and solidmechanics simulations, Technical Papers of SC2002, (2002).
[33] http://www.es.jamstec.go.jp