[共 同 研 究 成 果]

非構造格子に基づく流体数値計算法の並列化
中橋和博†、藤田健‡

1. はじめに「ＣＦＤが怒涛のごとく押し寄せてきた」、というのは２０年近く前に年輩の流体力学研究者から出た言葉である。確かに何百年ものあいだ紙と鉛筆、そして風洞実験で行われてきた流体力学研究にスーパーコンピュータという「飛び道具」を持って乗り込んで来た CFD （Computational Fluid Dynamics）は、今日では流体力学研究で非常に大きなウエイトを占めるまでになっている。その牽引役として航空宇宙分野も重要な役割を担ってきたが、すでに翼や機体形状あるいはエンジン要素の空力設計にＣＦＤは欠くことのできない道具になっている。そして、更に高い性能の航空機を開発するために、より高い精度でより早く結果をだす計算アルゴリズム開発の必要性が一層強くなっている。並列計算機はその要求に応えるための手段であり、かつ、長い歴史の流体力学研究に新たな展開をもたらしてくれるものでもある。
ここでは、CFD のこれまでの進展と近年の状況を簡単に述べるとともに、東北大学
情報シナジーセンターとの共同研究として実施した CFD ソルバーの並列化について、具体的なアルゴリズムとその結果を紹介する。
2. CFD の発展と現状
CFD のこれまでの進展を振り返った場合、いくつかのブレークスルーがあった。７０年代初頭に遷音速流れの計算が可能になり、７０年代中頃から８０年代にかけては境界適合座標の導入や陰的時間積分、高精度解法の開発等で翼型周りの流れの精度良い数値解析が飛躍的に進展した。今日の旅客機はマッハ 0.8 前後で飛んでいるが、このような音速近くを効率良く飛ぶための翼設計には CFD は不可欠なツールである。
９０年代の CFD は３次元実形状周りの流れの計算を目指した時期とも言えよう。
もちろん３次元計算は７０年代から行われている。しかし、多くは鈍頭円柱や単独の翼、あるいは翼と胴体程度の簡単な３次元形状の計算が主で、実際の飛行機周りの計算には数ヶ月、時には数年もの作業時間を要していたのも実状である。航空機の３次元形状設計にはとても使えるレベルではなく、CFD はきれいな絵を出すためのColorful Fluid Dynamics だとの悪口も言われていた。
流れの数値計算には、流れ場を非常に細かな網の目（格子）で覆い、個々の網の目毎に流速なり圧力を求める事を行う。その計算量は特に精度の要求される航空分野では膨大なものであり、７０年代から８０年代にかけては如何に効率良く計算を行うかが最重要な研究課題であった。そこで用いられたのが先に述べた境界適合座標の格子上での計算法である。網の目が順序よく並んでいるために何百万もの未知数をもった連立方程式でも効率よい計算が可能である。しかしながら、そのような整然と並んだ格子は幾何学的な制約を伴い、２次元では有効だったこの手法も、実際の飛行機全体†‡

東北大学大学院工学研究科航空宇宙工学専攻
東北大学大学院工学研究科航空宇宙工学専攻博士課程

形状の３次元計算では格子を作ることに何ヶ月もの時間を費やすこととなった。
９０年代は、それまで主流だった境界適合座標格子（構造的格子）に基づく計算法から離れ、実形状を扱うために格子点の並びに構造性のない非構造格子（三次元では４面体セルの集合）を用いる方法へとアルゴリズム研究者の関心が移った。しかしながら、構造性のある格子からの脱却には様々な新しい計算アルゴリズムの開発が要求される。３次元空間を何百万もの四面体で隙間無く自動的に埋めるにはどうするか、無秩序に並んだ格子点で流体運動の方程式を如何に効率良く解くか、あるいは航空機設計に必要な計算精度を非構造な格子で得るにはどうするか、等が主要な研究課題であった。
現在、この非構造格子を用いた計算法はかなりのレベルにまで達している。我々も非構造な格子生成法[1-3]や計算時間を大幅に短縮する陰的時間積分法[4]、物体移動を扱うための非構造格子オーバーセット法[5]等の独創的なアルゴリズムの提案を行ってきた。また、それら計算手法を用いて様々な実形状流れ問題に適用している[例えば 6-8]。
図１は、航空宇宙技術研究所で現在進められている超音速実験機プロジェクトに関する共同研究の一部である。これは次世代超音速旅客機開発の基礎研究として、CFDで低抵抗な形状を作りだし、それを飛行試験で実証しようというものである。抵抗の少ない形状を生み出すには、翼面上に設計者の希望通りの圧力分布を生み出す形状を求めるなり、あるいは抵抗値を最小化する形状の探索が必要となる。そのためにはCAD で定義された形状から素早く CFD 計算を行うことが求められる。図のような計算も、１０年前は数ヶ月を要していたが、数年前は数週間、そして今日では１日で数値解を得るまでに達した。並列計算機への対応により、更に数時間、数十分に収まるよう計算コードを改良することが今回の目的である。
3. 並 列 計 算 機 へ の 対 応
１万分の１のオーダーで抵抗係数を議論する航空機空力設計では、目に見えない微妙な形状変形でも空力性能を大きく変え、それが構造設計、ひいては航空機重量に影響を及ぼす。空力と構造等との厳しい妥協点探索のためにも膨大なシミュレーションが必要となり、数値計算時間の短縮は非常に重要である。
従来、CFD における計算アルゴリズムはベクトル計算機を念頭に開発されてきた。
そのプログラムの並列化は、しかしながら一様ではない。先に述べた境界適合格子（構造格子）上でのソルバーの並列化は比較的単純である。プログラム内の DO ループを並列処理して行くだけでかなりの並列化効果が得られる。一方、同様な方法で非構造格子上でのソルバーを並列化してもなかなか並列性能が上がらない。非構造格子上でのソルバーでは、無秩序に並んだ格子点の相対位置を覚えるために階層的な配列を多用しており、このことが並列処理の性能向上を妨げているようであるが詳しい原因は不明である。
そのため、ここでは計算領域を予め分割し、それぞれの領域を各 CPU に処理する
領域分割法に基づく並列化プログラムを作成した。領域間は Message Passing Interface（MPI）ライブラリを用いて情報交換を行う並列化手法である。このような領域分割に基づく並列処理は、そのための大幅なプログラム修正が必要であるが、一旦並列化プログラムができあがると、コンパイラに依存しない汎用性のあるコードともなる。

4. 領域分割法
領域分割法では、最初に計算前処理として計算空間に張られた格子の領域分割を行わなければならない。本研究では、ミネソタ大学で開発されたフリーソフトウェア
METIS[9,10]をベースに領域分割プログラムを開発した。図 2 に領域分割のフローチ
ャートを示す。最初に、METIS によりユーザーの設定した分割数で非構造空間格子を分割する。オリジナルの METIS では単一の要素からなる格子が基本となっているが、ここでは四面体やプリズム、ピラミッド要素が含まれるハイブリッド格子も取り扱えるように修正した。次に、領域間の接合面上の節点は複数の領域にまたがって存在するため、最も節点数の少ない領域に存在する節点を”Sending Vertices”と分類し、その節点に対応する他の領域の節点を”Receiving Vertices”とする。ここで”Sending Vertices”
は物理量や勾配などの情報を他の領域に送る点で、”Receiving Vertices”は他の領域からそれらの情報を受け取る点である。最後に”Sending Vertices”に分類された節点にオーバーラップする要素を追加し、その時追加された節点を”Receiving Vertices”に分類する。なお、その節点に対応する他の領域に存在する節点は”Sending Vertices”に分類される。
この方法により領域間の節点数および辺数の差が数％程度に抑えられ、CPU 間のロードバランスが整えられることが確認された。
5. 並列非構造格子 CFD ソルバー
並列化のベースとなるソルバーは、セル節点有限体積法に基づく非構造格子 3 次元
Navier-Stokes/Euler ソルバーを用い[4]、各領域の流体計算を別々の CPU に割り当てることで並列実行している。図 3 に情報交換の手順を示す。始めに、”Sending Vertices”にある情報（物理量や勾配）を単一の配列 sbuf に記録し別の PE に送る。sbuf が送られた PE ではその情報を配列 rbuf として受け取り、rbuf の情報を”Receiving Vertices”に割り振る。この方法は多数の節点の情報を 1 つの配列の通信で受け渡しするため、点から点へ情報交換を行うよりも通信のオーバーヘッドを減らす効果がある。なお、
このソルバーはベクトル化されているため、キャッシュベースの並列計算機のみならずベクトルベースの並列計算機についても利用することが可能である。
6. 計算精度とスケーラビリティ
計算精度を検証するため、研究室所有の ALPHA クラスターマシンを用いて ONERA
M6 翼に対する Euler 計算を行った。この空間格子の節点数は 234,621、 四面体数は1,210,089、辺数は 1,507,152 である。この空間格子を 16 分割した際の格子を図 4 に示す。図 5 に、ONERA M6 翼の 44%半スパン長における Cp 分布を、並列化されていな
い非構造 Euler ソルバーおよび実験値と比較する。並列化されていないものと完全に一致しており、並列化による計算精度の影響が無いことを確認した。
次に東北大学情報シナジーセンターのスーパーコンピュータ NEC SX-4 および汎用サーバーNEC TX7/AzusA を用いてスケーラビリティの測定を行った。前者はベクトルベースプロセッサを搭載し、後者はスカラーベースの並列計算機である。
始めに、航空宇宙技術研究所（NAL）が開発を行った超音速ロケット実験機（NEXST-1）に対する Euler 計算を、並列ソルバーを用いて行った。このモデルには

ピトー管やロケットブースターとの取り付け金具、発射台への固定金具など搭載され、実機に忠実な形状となっている（図 6 参照）。この形状は 3 次元 CAD ソフトにより定義され[8]、STL ファイルを出力し、表面格子生成プログラム Edge Editor[3]およびDelaunay 分割法[1]により空間を四面体で離散化した。このような複雑な形状に対しても非構造格子を用いることで容易に格子生成を行える。この非構造空間格子の節点数は 1,016,066、辺数は 6,715,503、四面体要素数は 5,524,311、三角形要素数は 350,254である。さらにこの空間格子を 7 つの領域に分割したものを図 7 に示す。この Euler
計算に必要なメモリサイズは約 2.2GB であり、NEC SX-4 を用いて並列化性能を測定した（図 8 参照）。十分な並列化性能を達成しているものの、20PE からわずかに頭打ちになっている。全メモリサイズは 1PE の計算と大きく変わらず、32PE の時に 1PEあたりのメモリ使用量が 70MB 程度になってしまうため、十分なベクトル長が確保できず性能を十分発揮していないことが原因である。実際、図 9 に 1PE あたりの平均ベクトル化率と MFLOPS 値を示すが、32PE の際に平均ベクトル化率が 98%を下回っている。なお、図 9 のエラーバーは各値の最小値および最大値を示す。
次に、より大規模なスケーラビリティの測定のために、図 10 に示すリージョナルジェット機に対する Navier-Stokes（NS）計算を行った。この格子では、境界層を十分離散化するため壁面付近をプリズム要素で離散化するハイブリッド格子を用いており、さらに主翼後縁部分に存在するわずかな厚みも再現しているため非常に大きなメモリサイズを要する。空間格子の節点数は 2,180,582、辺要素数は 10,577,008、四面体要素数は 3,839,284、プリズム要素数は 2,943,184、ピラミッド要素数は 27,322 となっており、計算に必要なメモリサイズは約 4GB である。
NEC SX-4 および TX7/AzusA におけるスケーラビリティの結果をそれぞれ図 11、図12 に示す。流体ソルバーは陰解法を用いているため 3000 ステップ程度で収束解を得ることができるが、この 3000 ステップの計算を SX-4 の 1PE で実行すると 3 日かかる。
一方、32PE で実行することでわずかに 2 時間半で結果を得ることができる。これらにより、並列非構造格子ソルバーがベクトルベースおよびスカラーベースの並列計算機で優れた並列化性能を達成することが示された。なお、SX-4 の 32PE を用いることで 10GFLOPS もの性能を達成した。
7. おわりにMETIS をベースとした領域分割プログラムと、MPI ライブラリにより並列化された非構造 CFD ソルバーを用いた並列非構造格子 CFD 計算システムを構築した。このソルバーを情報シナジーセンター所有の NEC SX-4 および TX7/AzusA を利用し、並列化性能の測定を行った。その結果、ベクトルベース、スカラーベース双方の並列計算機で優れた並列化性能を達成し、計算精度も十分であることが示された。特にベクトルベースの NEC SX-4 で得られた高い並列化性能は重要である。
非構造格子を用いた流体計算法は実形状を容易に扱えることが大きな利点であり、例えば図 13 のような、表面に 400 程ものディンプルを持ったゴルフボールの空力解析や[11]、図 14 のような昆虫を忠実に計算機内で再現し[12]、その羽ばたき飛翔を調べるような計算も可能になりつつある。しかし、リアルな形状の再現には膨大な計算量も必要である。今回の並列化ソルバーの開発により、今まではとうてい無理であった様々な流体問題にチャレンジして行くことが可能となる。航空機開発においても、

これまでの静的な解析だけでなく、飛行の動的な解析までをも視野に入れた数値シミュレーションが始まろうとしている。今後の超大規模な流体数値計算が楽しみである。

8. 謝辞
本研究は、東北大学情報シナジーセンターとの共同研究として行われた。SX-4 お
よび TX7/AzusA 上での計算コードの並列化・最適化および実行に際して情報シナジーセンターから適切かつ有益な助言を頂いた。この場を借りて感謝したい。
9. 参考文献[1] Sharov, D. and Nakahashi, K., “A Boundary Recovery Algorithm for Delaunay
Tetrahedral Meshing,” 5th International Conference on Numerical Grid Generation in
Computational Field Simulation, pp.229-238, 1996.
[2] Sharov, D. and Nakahashi, K., “Hybrid Prismatic/Tetrahedral Grid Generation for
Viscous Flow Applications,” AIAA Journal, Vol.36, No.2, pp.157-162, 1998.
[3] Ito, Y. and Nakahashi, K., “Direct Surface Triangulation Using Stereolithography Data,”
AIAA Journal, Vol.40, No.2, pp.490-496, 2002.
[4] Sharov, D. and Nakahashi, K., “Reordering of Hybrid Unstructured Grids for
Lower-Upper Symmetric Gauss-Seidel Computations,” AIAA Journal, Vol.36, No.3, pp.
484-486, 1998.
[5] Nakahashi, K., Togashi, F., Sharov, D., “An Intergrid-Boundary Definition Method forOverset Unstructured Grid Approach,” AIAA Journal, Vol.38, No.11, pp.2077-2084, 2000.
[6] Nakahashi, K., “Progress in Unstructured- Grid CFD”, Proc. 1st ICCFD, Computational
Fluid Dynamics 2000, Springer, pp.3-13, 2000.
[7] Murayama, M., Nakahashi, K., Sawada, K., “Simulation of Vortex Breakdown Using
Adaptive Grid Refinement with Vortex-Center Identification,” AIAA Journal, Vol.39,No.7, 1305-1312, July 2001.
[8] Fujita T., Ito Y., Nakahashi K. and Iwamiya T., “Computational Fluid Dynamics
Evaluation of National Aerospace Laboratory Experimental Supersonic Airplane inAscent,” Journal of Aircraft, Vol. 39, No. 2, pp. 359-364, 2002.
[9] G. Karypis and V. Kummar, “A fast and high quality multilevel scheme for partitioning
irregular graphs”, Technical Report 95-035, University of Minnesota, 1995. A shortversion appears in Intl. Conf on Parallel Processing 1995.
[10] http://www-users.cs.umn.edu/~karypis/[11] Kato T., Fujita T., Ito Y., Nakahashi K. and Kohama Y., “Computational Analysis ofSeparated Flow around a Golf Ball Using Unstructured Grid CFD," AIAA Paper
2001-2569, 15th AIAA Computational Fluid Dynamics Conference, 2001.
[12] Togashi, F., Ito, Y., Murayama, M., Nakahashi, K. and Kato, T., “Flow Simulation ofFlapping Wings of an Insect Using Overset Unstructured Grid,” AIAA 2001-2619, 15thAIAA Computational Fluid Dynamics Conference, 2, 2001.

Fig.1 Computed pressure distribution around the NAL experimental supersonic airplaneNEXST-2.

1. Division by MeTiSDivisionline

2. Decision of calculating region of vertices on boundary

3. Mesh overlapping

Receiving Vertices

Sending Vertices

Fig. 2 Flowchart of division of unstructured grid

PE1Sub-domain 1

Sub-domain 1

PE1

sbuf

rbuf

PE2

PE2
rbuf
Sub-domain 2

Sub-domain 2

sbuf

Fig. 3 Send/Receive data between two sub-domains

Fig. 4 Partitioned grid of ONERA-M6 wingExperiment-1.5

Parallel-1

No Parallel

Cp-0.500.5
10Chord Length1

Fig. 5 Comparison of Cp distribution at 44% semi-span

Fig.6 NAL NEXST-1 precise model

Fig.7 Partitioned grid of NAL NEXST-1

35Ideal
30Parallel Unstructured Solver
SPEEDUP25201510500816PE

2432

Fig. 8 Scalability result 1 on NEC SX-4

Ave. V. Op. Ratio (%)

100
99
98
9796
95
12572032

72032

PE500

Ave. MFLOPS

400
300
200
1000125PE

Fig. 9 Average vector operation ratio andaverage MFLOPS par PE

Fig. 10 Regional jet airplane35 Ideal
30 Parallel Unstructured SolverSPEEDUP

25201510500816PE

2432

Fig. 11 Scalability result 2 on NEC SX-4

16 Ideal
14 Parallel Unstructured SolverSPEEDUP1210864
20048PE

1216Fig. 12 Scalability result 2 on NEC TX-7/AzusA

Fig.13 Flow simulation around a golf ball.

Fig.14 Flow simulation around a hornet.