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

＊ここに掲載するのは、エジンバラ大学のAdrian Jackson博士およびSergio Campobasso博士によるHECToRレポート「OPTIMISING THE PARALLELISATION OF A HARMONIC BALANCE NAVIER-STOKES SOLVER FOR THE ULTRA-RAPID ANALYSIS OF WIND TURBINE, TURBOMACHINERY AND AIRCRAFT WING PERIODIC FLOWS, Adrian Jackson, M. Sergio Campobasso, University of Edinburgh, 29 March 2013」を要約したものです。

概要

COSAコードはFORTRANで記述され、定常、時間領域（TD）[1-3]およびハーモニックバランス（HB）[2]による流体ソルバーを含みます。定常およびHB方程式の積分は、マルチステージのRunge-Kutta $RK$ （スムーサー）を用いた陽的なマルチグリッド法をベースにしています。TDでの非定常問題の解法は、Jamesonの二重時間刻みステップ法を用いています。乱流効果はMenterの2方程式K-ωずれ応力輸送（SST）乱流モデル[1]を用いて組み込まれています。

COSAはLow-Speedプリコンディショナ[2,3]も用いており、水平軸風力タービン（HAWTs）の様な低速流にも利用可能です。3つのソルバー（定常、TB、HB）は、その計算量が膨大な計算リソースを必要とするため並列化されています。既存のコードは、OpenMP、MPI、Hybrid $mixed OpenMP and MPI$の3種の並列手法を用いています。

OpenMP並列は全てのソルバーで利用でき、マルチグロックグリッドでのブロックや、HBソルバーのハーモニクス、および少数のブロックやハーモニクスを使用する問題でのグリッドポイント上で異なる並列処理が可能です。

MPI並列はマルチブロックグリッドのブロックを分散し、シミュレーションをmぴプロセス間で分散処理します。ブロック間では通信が生じ、ブロックの端（COSAではcutsと呼ばれます）のデータは隣接ブロックとの通信（ハロ通信）を必要とします。MPI並列プロセス最大数はジオメトリーパーティション（グリッドブロック）数で制限されます。

ハイブリッド並列ではMPIは、シミュレーションの種類によってハーモニクスOpenMP並列、あるいはグリッドポイントOpenMP並列と組み合わされます。

COSAの既存の並列処理は、その効率が最適化されていません。OpenMP並列はデスクトップマシンの様なローカルコンピュータや小規模共有メモリークラスタでの並列処理は可能ですが、大規模リソースの利用は困難です。更にOpenMP並列は少数コアにおいても十分スケールするとは言い難い手法です。MPIは一般的に有効で改善可能なエリアが有りますが、利用可能なコア数がシミュレーションで用いるグリッドブロック数の最大値に制限されます。ハイブリッド並列はこの制限に対処するもので、ブロック数より多くのコア数を用いる事を可能にしますが、MPIのみの並列処理よりも効率がよく有りません。

よって以下の2つの部分についてHECToRで効率を上げることを目標にしました：

• HBソルバーのMPI通信の最適化


• HBソルバーのOpenMPおよびハイブリッド並列の最適化

既存コードの性能

MPIコードはOpenMPコードよりも優れていますが、特に大きなプロセス数で最適化が十分では有りません。MPIコードはシリアルコード機能に一対一に対応づけられて設計されており、特定のグリッドcut（隣接ブロック間の境界）に関連する全てのハロセルデータを一度に送るのでなく、一度に一つのハロセルデータを転送しています。

効果的なMPI実装では、ブロック総数を変化させずにブロックサイズを増やすと、MPI効率に対して無視出来ない影響を及ぼします。何故なら、MPI通信の回数は殆ど変化せず、計算と通信に掛かる時間比はほぼ一定になるからです（これはMPIメッセージサイズと同時に各ブロックにおける計算量が増加するためです）。

逆にこのMPIコードの中では、ブロックサイズの増大はMPI並列効率を劣化させます。何故なら、この操作は小さなMPIメッセージを大量に発生させるからです。理想的な並列コードではメッセージ数は一定になり、そのサイズはメッセージ数に依存して増えるのでなく、ブロックサイズに依存します。以下の表に例を示します。これはHECToRフェーズ2b上での、既存のMPI並列HBソルバーのスケーリングテスト結果です。テストケース2を用いてベンチマーク実行しました。予想通りに、大量の小さなMPIメッセージが発生する細密グリッドの場合より、粗いグリッドの方がMPI並列効率は良好です。この問題を除けば、既存のMPI HBコードは大規模問題でも十分な性能を持ちます。

COSAの実行時間の約80％を占めるのはユーザルーチン（COSAのサブルーチン）で、MPI通信は15％を占めています。また、10個のサブルーチンがコード全体の約60％を占めています。よってこれらのルーチンをOpenMP性能とMPI通信部分の最適化対象とします。

MPI並列の最適化

既存のMPI実装では、ハロデータ通信あるいは隣接プロセスとのcutデータ通信において、膨大なメッセージを送受信します。これは既存コードが一度に、隣接プロセスへ境界データの小さな一部をおくるためです。

各要素に対してメッセージを送信せず、送信データあるいは受信したデータを配列へ収集して送信するように改変しました。

結果として、最適化されたMPI通信手法がコード性能を改善していませんでしたが、新コードは既存コードに比べ、メッセージ送信回数が2桁少ないことが示されました。

ここで、ノード当たり 2 x 8 コアの Intel Xeonプロセッサ、および Infiniband QDR interconnectで構成された別のスーパーコンピュータでは、既存のコードよりも2桁高速化されていました。
結論としては、HECToRネットワークは小さなメッセージに対して良好な性能を持ちますが、ここで実装されたメッセージ集約はHECToRでは改善されませんでした。しかしながら他のシステムによっては性能改善されます。

集団通信の最適化

MPI集団通信のメインルーチンであるcombineforcesは、ハーモニクス数のループ内から呼ばれており、all reduceがnbody*$(2*nharms$+1)回呼ばれています。

そこで、combineforcesが呼ばれる場所を移動して、all reduceで通信されるデータを集約して、全てのall reduceをcombineforcesの呼び出し毎に一回のall reduceへ集約したところ実行時間が1桁小さくなりました。

MPI-I/Oの最適化

COSAはシミュレーションの最後に1回リスタートファイルを、同じくハーモニクス当たり1つのflowtecファイルを書き出します。大規模シミュレーションではリスタートファイルは数GB、flowtecはGBサイズになります。既存コードでは、MPI I/Oが装備されていますが、個別のデータ要素に対して実行しています。

MPI I/Oは一回の呼び出しで大量データをファイルに書き込む際に性能を引き出すことができるものであるため、書き込むデータを配列へ集約して一度に書き込む様に修正しました。

リスタートファイルが350MB、31個のflowtecファイルが各々20MBのテストケースの場合に測定を行ったところ、全実行時間はI/Oを含めると既存コードはI/O無に比べて60％増加していましたが、新コードでは30％の増加に削減されました。

これ迄の最適化を全て含めたMPIコードは、テストケース1の場合に256と512コアで実行時間は約半分になりました。

行列ベクトル積の最適化

COSAは、陽的なfully explicit Runge-Kutta（FERK）スムーサーと陰的なPoint-Implicit Runge-Kutta（PIRK）スムーサーの2つのモードを搭載しています。PIRKではHBソース項は各Runge-Kuttaステージにおいて陰的に解かれます。こうして各ステージで、次元が$npde*((2*nharms$+1))のncell個の線形システム解が要求されます。PIRKによるHB時間積分は、通常のFERKのHB時間積分に対して解の時間マーチングに大きな擬似時間ステップが数値安定性のために用いられます。PIRKは実際にはMGサイクルの繰り返し数がFERKより遥かに少なく収束します。

このレポー内トの全ての時間計測は比較的時間ステップの少ないテストケースのためFERKを用いていますが、現在COSAの実行で頻繁に用いられるのはMGサイクル数がより少ないために高速なPIRKです。そこでLAPACKによるPIRKの性能を検証しました。その結果、512プロセスで、LAPACKを用いた新バージョンの実行時間は既存コードに比べ約半分になりました。

ハイブリッド並列処理の最適化

既存のハイブリッドコードは、シミュレーションタイプによって、ブロックやハーモニクスあるいはグリッドセル上で並列化されています。ここではハーモニックバランスに着目して最適化した結果を報告します。

ハーモニックバランス部分のOpenMP並列領域はそれほど多くの実行時間はかかっていませんでしたが、ファーストタッチによる初期化や、OpenMP呼び出し回数が多くオーバーヘッドが大きいため、余分な並列領域の削除や外側ループへのスレッド並列の移動などの再設計をしました。

次にOpenMP並列化されていない、I/OとMPI通信を最適化しました。

MPI通信はcutとハロおよびデータのループに対して実行されます。ここでcutは独立なのでスレッド並列が可能です。しかしながらこれらはMPI通信を含むため、ここでは通常のMPI_INITではなく、スレッド版MPIライブラリのMPI_INIT_THREADを用います。またここで使用するMPIライブラリは、個々のOpenMPスレッドがMPI通信が可能なMPI_THREAD_MULTIPLEをサポートするものです。

/OはMPI I/Oを用います。一般に、I/O操作はブロック間およびブロック内のハーモニクスに関して独立ですが、ファイルオープンとクローズで全プロセスの操作が含まれる集団通信が存在します。そこで、OpenMPスレッドがリスタートファイルやflowtecファイルに独立に書き込める様に修正しました。最終的な結果として、MPIタスク当たり1ブロック、4スレッドで、ハイブリッドコードは4倍にコア数を増やすと性能は3.68倍になりました。

まとめ

	MPI$512コア$	旧ハイブリッド$2048コア$	新ハイブリッド$2048コア$
実行時間[秒]	1693	772.16	520
加速率	-	2.19	3.26

表:MPI性能を基にした2048コアを用いたハイブリッド性能の新/旧コード比較

最終結果として、既存コードが4倍の計算リソースで半分の実行時間、つまり効率が約55％だったのに対し、改善後の並列効率は90％になりました。つまり同じ問題を半分の時間で実行出来る様になりました。

謝辞

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

文献

[1]	M.S. Campobasso, A. Piskopakis, M. Yan, Analysis of an Oscillating Wing in a Power-Extraction Regime Based on the Compressible Reynolds Averaged Navier-Stokes Equations and the K-ω SST Turbulence Model, ASME paper GT2013-94531, presented at the ASME/IGTI Turbo Expo 2013 Technical Conference, 3rd-7th June 2013, San Antonio, Texas, USA..
[2]	M.S. Campobasso and M.H. Baba-Ahmadi, Analysis of Unsteady Flows Past Horizonatal Axis Wind Turbine Airfoils Based on Harmonic Balance Compressible Navier-Stokes Equations with Low-Speed Preconditioning, ASME Journal of Turbomachinery, Vol. 134, no. 6, 2012. DOU: 10.1115/1.4006293.
[3]	M.S. Campobasso and J. Drofelnik, Compressible Navier-Stokes analysis of an oscillating wing in a power-extraction regime using efficient low-speed preconditioning, Computers and Fluids, Vol. 67, 2012, pp. 26-40. DOI:10.1016/j.compfluid.2012.07.002.
[4]	http://www.scalasca.org/
[5]	A. Jackson and M. S. Campobasso. Shared-memory, distributed-memory, and mixed-mode parallelisation of a cfd simulation code. Computer Science - R&D, 26$3-4$:187.195, 2011.
[6]	A. Jackson, M. S. Campobasso, and M. H. Baba-Ahmadi. On the Parallelization of a Harmonic Balance Compressible Navier-Stokes Solver for Wind Turbine Aerodynamics, pages 747.761. ASME International, 2011.