プログラムの高速化・並列化サービスの事例

＊ここに掲載するのは、STFC Daresbury LaboratoryのXiaohu Guo博士らによるHECToRレポート「Massively Parallel Computing for Incompressible Smoothed Particle Hydrodynamics, Xiaohu Guo$a$, Benedict D. Rogers$b$, Peter. K. Stansby$b$, Mike Ashworth$a$, a)STFC Daresbury Laboratory, b)University of Manchester, March 2014」を要約したものです。

概要

射影ベース圧力補正付きISPH（incompressible smoothed particle hydrodynamics）法が内部流れに対して高い正確性と安定性を持つことが示されました。
このプロジェクトでは、沖合と沿岸構造における激しい自由表面流に対するISPH法を用いた流体ソルバーの最適化を行いました。ISPHのこれまでのベンチマークから、最近傍探索と圧力ポアソン方程式ソルバー（PPE)が主要なコストを占めることが示されています。並列化においては、ヒルベルト空間充填曲線（HSFC)とZoltan[5]パッケージが領域分割に用いられ、その空間的局所性が最近傍探索の性能を決めます。

主要な作業は以下の通りです：

・粒子とその物理場データブロックを適切なパーティションへ送る機能を提供するマップカーネルを修正し、MPI片方向通信とZoltanの分散ディレクトリユーティリティで最適化しました。これによりマップカーネルの占有割合は、粒子が領域を跨ぐ様に均一に分散されるシミュレーションで20％以下に削減されました。

・最近傍探索カーネルは、前処理付き動的ベクトルとリンクリストを用いて書き直されました。さらに以下の最適化を行いました：
-　ジェネリックインターフェイスを追加して入力パラメータ設定から呼び出しを可能にしました。
-　粒子間距離を保存せずに常に計算するようにしたところ、キャッシュ効率が向上し、メモリーが縮小されました。これにより粒子数が小さい場合は10倍、粒子数が大きい場合は3倍まで高速化しました。
-　平滑化カーネル関数の最適化によりメモリーが削減されました。

・HSFCでのセルの再順序化が実装されましたが性能改善には至りませんでした。原因は再順序化処理そのものでした。これによる全体性能の改善は5％以下でした。

・圧力ポアソン方程式（PPE)ソルバーに対して以下の最適化を行いました：
-　new Yaleスパース行列フォーマットをPETScのCSRフォーマットへ置き換えました。これにより大幅にメモリーが削減されました。
-　全体行列構築前に粒子を再番号付けました。その結果、全体行列への値の挿入操作が局所操作になり、PPEソルバーの性能は大きく向上しました。PETSｃのマルチグリッド・プレコンディショナを用いると、コードは100M個の粒子に対して8千コアまでスケールします。

プロジェクトの概要

このISPHコード[1,2]は沖合および沿岸における工学的な課題を解くために数年に渡って開発されてきました。それは例えば、複雑で高度に非線形性を持ち無秩序な波の運動であり、そこでは砕波、海嘯の伝搬、空気混入、構造相互作用、衝撃などが含まれます。当初のコードはシリアル版でしたが、2010年8月以来、Manchester SPH groupとDaresbury Application Performance Engineering Group$APEG$のEPSRCジョイントプロジェクトにおいて開発が進められました。コードはFortran90で書き換えられMPIで並列化されました。ここではZoltanを用いた領域分割と動的負荷バランスおよびPETSｃのスパース線形ソルバーが用いられました。

作業項目1：粒子マッピングカーネルの最適化

ISPHは領域分割にZoltan[http://www.cs.sandia.gov/zoltan/Zoltan.html]を利用します。これまでは粒子をセルに配置してから、それを用いて隣接リストを作成していました。現在は全体セルの分割にヒルベルト空間充填曲線を用いた結果としてデータをプロセッサーへ割り当てます。このマッピングルーチンは粒子とその物理場データブロックを適切なパーティションへ送る機能を提供します。1024MPIタスクを用いたCrayPAT解析から、主要なボトルネックはマッピング関数に含まれるMPI_GATHERVにあり、全MPI通信の86％を占めていることが判りました。これは、このマッピング関数が、個々のセルの識別子：cell_idとパーティション識別番号：IDをグローバルに保持していることが主たる原因です。ダムの崩壊のような高度に乱雑な流れの場合には、通常はセル数が粒子数よりも遥かに大きくなります。粒子数が大きい場合は、グローバルに共有されるこのマップの入れ替えは非効率になります。この問題について以下の2種の方法を提案しました：

·回避策1：粒子マッピングカーネルのMPI片方向通信による最適化

各時間ステップではマッピングルーチンを一回だけ更新するようにします。MPI_PUTのMPI派生タイプを用いれば、MPIタスク当たりのMPI_PUT操作数を削減可能になります。

·回避策2：Zoltan分散ディレクトリの利用

このようなセル数が粒子数よりも何倍も大きい場合は、こうしたグローバルマップの維持と更新は極めて高コストです。ここで、Zoltanの分散ディレクトリ機能を用いることを検討しました。これは、データマイグレーション後のセル配置管理やグローバルにアクセス可能な情報の作成に対してランデブーアルゴリズム$rendezvous algorithm$[3]を用います。分散ディレクトリは、メモリーとプロセッサ時間の意味での負荷を均衡させてマッピングのボトルネックを解消します。

残念ながらZoltanライブラリは分散ディレクトリにFrotranインターフェイスを装備していないため、ここでFortranインターフェイスを作成し、これを、2048MPIタスクを用いた55M粒子のダム崩壊テストケースに適用しました。
全体として20%以上のコスト削減が示されました。

作業項目2：隣接リスト検索モジュールの最適化

一般に、スピードや並列化に関わらず、大規模粒子数のシミュレーションは効率的な隣接探索アルゴリズムが実装されている場合のみ可能であることは明らかです。通常は隣接リスト作成には2つの方法がとられます：第一の方法はセルリンクリスト法（CLL)です。これは、セルとその隣接セル内に置かれた粒子間の相互作用を計算するのに用いられるすべてのセルへのリンクを示すリストを作成します。第二の方法はベルレリスト（Verlet List $VL$）と呼ばれるもので、各粒子へのリンクリストを作成します。ISPHコードはVL法よりもメモリーが少なくて済むCLLを用います。リスト作成の最初のステップでは、粒子はセル内に配置されるため両方の方法で同じです。リンクリスト手法のような隣接リストを作成するには、粒子をセルに配置するやり方が複数ありますが、前処理付き動的ベクトル（preconditioned dynamic vector）[4]に対してはどれも似た性能とメモリー効率を持ちます。そこで、前処理付き動的ベクトルとリンクリスト法で隣接リストモジュールを書き直しました。リンクリスト法を用いれば、アダプティブ粒子シミュレーション（粒子にはアダプティブメッシュリファインメント（AMR）と同様に分割と融合が行われます）に適切なデータ構造を持つことが可能になります（これは、リンクリスト法を用いない場合には各粒子の隣接リスト用のメモリの事前確保が困難になるためです）。前処理付き動的ベクトを用いた隣接検索モジュールは、中規模粒子数の非アダプティブSPHシミュレーションに用いることが出来ます。

既存コードでは、各時間ステップで粒子間距離は一度だけ計算され、3次元配列に格納されます。各配列サイズは粒子数×リンク数で、繰り返し多くの浮動小数点演算を要します。新しいコードでは、この配列をメモリーがより少なく、浮動小数点演算がより大きな関数に置き換えます。これにより少数粒子では10倍、大規模粒子では3倍の高速化が示されました。

通常は、粒子の番号は最初に与えられた順序で番号付けされています。このメモリーレイアウトはシミュレーション中は固定です。粒子の再番号付けにより、メモリーアクセスパターンが変更され、大規模データ内の小さなサブセットをキャッシュに確保することが可能になります。このキャッシュ内のデータブロック上で計算を行うようにします。キャッシュ内のデータを使い回して計算が可能であれば、メインメモリーへのアクセスを削減できます。
SPHの隣接検索法は、セルの並びに従って演算するので、セルの順序に粒子を並べる方が自然です。ISPHコードでは、セルは空間充填曲線上に並ぶので、粒子もまたその順序に従って並べられます。結果として、こうした再順序化と粒子の移動に伴って、メモリーアクセスパタンは数時間ステップ毎に変化します。また、これによりHSFCはデータの局所性を維持することから、キャッシュ効率が改善します。

作業項目3：圧力ポアソンソルバーの最適化

既存のシリアル版ISPHコードにおいて、全計算時間の約47%が圧力ポアソンソルバーで占められています。その効率化のためにPETSｃ[6]ソフトウェアをISPHコードへ組み込みました。
シリアル版のISPHはスパース行列にnew Yaleフォーマットを用いますが、デフォルトのPETSｃの行列表現には一般スパースAIJフォーマット（Yaleフォーマットやcompressed sparse row format（CSR)とも呼ばれます）を用います。そこで解くべき行列のPETSｃ行列への変更が必要です。デフォルトのAIJフォーマットの内部行列表現はゼロベースインデックスに従います。

PETSｃは行列を連続した行で分割しますが、ISPH内では各行は特定の粒子に隣接する全ての粒子を意味します。再番号付けされた行列が用いられるため、各パーティションはそれ自身の行列を構成することになります。HSFCが用いられる場合、解くべき線形系は元のものよりも狭いバンド幅を持つため、PPSソルバーの性能が向上します。
既存のISPHコードでは、ヤコビ法プレコンディショナと安定化BiCGS法を用います。PETScによる改善後は、HYPRE BoomerAMGやPETSc GAMGなどのマルチグリッド・プレコンディショナを用いることで収束性やロバスト性および実行時間の改善が示されました。

·ベンチマークと性能解析

テストケースはダム崩壊と静水です。全粒子数を2Mから100Mまで変化させます。HECToRのCrayXE6システム上で実行します。これは2816ノードで、各ノードは2つのAMD 2.3GHz　16コアプロセッサーから構成され、全体で90,112コアです。理論ピーク性能は800Tflops強です。ノード当たり32GBのメインメモリーを搭載し、32コアで共有し、全体で90TBです。プロセッサはCray Gemini通信チップにより高バンド幅の3次元トーラスで結合されています。

HSFCを用いたZoltan領域分割の時間は、1024コア上では1.01%です。コア数がより小さい場合はさらに小さくなります。

HYPRE BoomerAMGプレコンディショナとGMRESを用いた場合のISPHの並列効率を、1024コアから8192コアまで計測しました。行列構築時間は実際のソルバーに掛かる計算時間に比べ極めて小さいです。PPEソルバーの平均並列効率は85%以上で、8192コアまでスケールしています。

MPIはISPH全体に対して依然として41%を占めています。主要なコストはZoltanに起因するMPI_Allreduceです。ダム崩壊ケースでは空のセルが存在しており、将来これを除けば、セル数が削減されてZoltanの性能は向上すると見込んでいます。PPEソルバーの割合は43.4%となり、主要なコストはHYPRE BoomerAMGです。カーネル計算は隣接リストの最適化により2%以下になりました。

粒子数を2Kから14Kまで振って、粒子の再番号付けの性能を調べました。粒子数に依存しますが、4.4%から8.5%の性能改善が示されました。性能向上が期待ほどでない理由には、再番号付け処理自身のオーバーヘッドがあります。

この仕事の一部は、8th International SPHERIC SPH Workshop in Trondheim, June 2013 にて報告されました。

謝辞

このプロジェクトは、nAG Ltd.が運営するHECToRの分散計算科学および工学$CSE$サービスの基に実行されました。英国の国立スーパーコンピューティング・サービスである、HECToR：英国リサーチ・カウンシル・ハイエンド計算サービスは、リサーチ・カウンシルを代行するEPSRCが管理しています。そのミッションは英国学術界の科学および工学の研究支援です。HECToRスーパーコンピューターは、UoE HPCx Ltd.およびnAG Ltd.のCSEサポートサービスにより管理運営されています。

文献

[1]	Lind, S.J., Xu, R., Stansby, P.K. and Rogers, B.D. Incompressible Smoothed Particle Hydrodynamics for free surface flows: A generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves, J. Comp. Phys, 231:1499-1523, 2012
[2]	Guo X, Lind SJ, Rogers B D, Stansby P K, Ashworth M. Efficient Massive Parallelisation for Incompressible Smoothed Particle Hydrodynamics with 108 Particles. Proc. 8th International SPHERIC Workshop, Editor J.E. Olsen. pp 397-402. 4-6 June 2013
[3]	A. Pinar and B. Hendrickson, Communication Support for Adaptive Computation, in Proc. SIAM Parallel Processing 2001
[4]	J.M. Dominguez et al. Neighbour lists in smoothed particle hydrodynamics, Int. J. Numer. Meth. Fluids, 2010. DOI: 10.1002/d.2481
[5]	Erik Boman, Karen Devine, Lee Ann Fisk, Robert Heaphy, Bruce Hendrickson, Courtenay Vaughan, Umit Catalyurek, Doruk Bozdag, William Mitchell, James Teresco, Zoltan 3.0: Parallel Partitioning, Load-balancing, and Data Management Services; Developer's Guide, Sandia National Laboratories, 2007, Tech. Report SAND2007-4749W
[6]	PETSc Manual Page, http://www.mcs.anl.gov/petsc/documentation/index.html