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

＊ここに掲載するのは、エジンバラ大学EPCCのSavvas Petrou博士によるHECToRレポート「dCSE Fluidity-ICOM: High Performance Computing Driven Software Development for Next-Generation Modelling of the Worlds Oceans, Xiaohu Guo1), Mike Ashworth1), Gerard Gorman2), Stephan Kramer2), Matthew Piggott2), 1)Daresbury Laboratory, 2)Imperial College London」を要約したものです。

概要

Fluidity-ICOM $Imperial College Ocean Model$プロジェクトは、2009年7月から2010年5月まで行われました。Fluidity-ICOMはオープンソースです（https://launchpad.net/fluidity）。
このプロジェクトは、プロファイリングと最適化、メッシュアダプティビティのスケール化、並列I/Oの３部から成ります。全体目標は、大規模スーパーコンピュータ上においてFluidity-ICOMを効率的に動作させることです。

Fluidity-ICOMは、アダプティブ非構造有限要素CFDコードFluidityをベースに構築されています。これは3次元非静水圧並列マルチスケール海洋モデルで構成され、非構造異方性アダプティブメッシュ上に様々な有限要素法と有限体積法を実装しています。よってネストされたグリッドを必要とせずに、非常に広い範囲の結合された構造を単一の数値シミュレーションで正確かつ効率的に表現することができます$Ford et al., 2004; Pain et al., 2005; Gorman et al., 2006$。

非構造メッシュを用いた海洋モデル開発は、これまでの有限差分モデルに比べ大きく複雑性が増します。モデルの数値計算コアは別として、これまで標準的なものがないプリポストツールの開発に多くの工数が必要となります。例えば、初期のメッシュ生成は複雑な海底地形や海岸線を表現できる必要があり、これは海洋のダイナミクスに重要な影響を及ぼします。 実際には、離散化された領域がどれほど現実を表すのかと、限られた計算資源で数値シミュレーションを行うのにどれほど適しているかのトレードオフがあります。これを実現するために、CADモデルを使用する標準パッケージには適さないジオメトリを生成する、Terreno（Gorman et al.,2008）などの特殊なメッシュ生成ツールが開発されています。

Fluidity-ICOMはまた、異方性メッシュアダプティビティを用いたモデリング誤差推定を制御する数値メッシュ最適化が可能です。メッシュアダプティビティにより大きな負荷分散が予想されることから、動的不可バランス手法が必須です。メッシュの再バランス化を行う場合、拡散再分割法と領域再マッピングによる並列再分割法の組み合わせを使用してデータの移動を最小限に抑えます。 これは現在、ParMETISとデータ移動のための特注コード（Gorman et al.,2009）の組み合わせを使用して行われていますが、代わりにZoltan（http://www.cs.sandia.gov/Zoltan/）を組み込む作業が行われています。

スパース線形方程式系の求解はFluidity-ICOMにおいて最も時間が掛る部分です。Fluidity-ICOMは、条件数が少ない場合はPETSc（http://www-unix.mcs.anl.gov/petsc/petsc-2/index.html）による標準的なプレコンディショナとソルバーを用います。条件数が多い場合はAMGプレコンディショナを用います。これは海洋問題に特化して設計されたもので、通常のプレコンディショナやブラックボックスAMGプレコンディション試験されたものを上回る性能を示します。しかしながら、AMGの標準的なプロセスは粗く、スケーリングが難しくなるため常に注意を払う必要があります。

Fluidity-ICOMは利用可能である場合は、最新の標準的な外部ソフトウェアライブラリを利用することが可能です。例えばPETScはスパース線形系の求解に、またZoltanは様々な並列データ操作に用いられます（これらはともに互換性のあるオープンソースライセンスを持ちます）。Fluidity-ICOM実行時のユーザ定義関数や診断ツール、問題設定に対してはPythonが用いられます。総数17の外部ソフトウェアパッケージと3種の言語（Fortran, C++, Python）が用いられています。

ベンチマークデータ

最初に、初期実行での不安定性を落ち着かせてからチェックポイントで保存し、その後、ベンチマーク実行を行うようにしました。
ベンチマークによっては、上述の最初の初期実行時間が非常に長くなることがあります（たとえば、1年間の還流シミュレーション時間など）。そこで、初期実行では中程度の64プロセスで実行します。各プロセスが自身のパーティションを読み書きするため、結果として得られるチェックポイントは予め64のサブ領域に分解されていることになります。次に、Fluidity-ICOMツールキットプログラムであるfreedecompを使用して、様々な数のサブ領域に問題を再分割し、強スケーリング解析を実行します。通常、強スケーリング解析では、パーティションの数の範囲は問題全体のサイズによって制限されます。与えられた問題に対してプロセッサが多すぎると、ハロのサイズがサブドメインのサイズに近づいていくため、計算コストに比べて通信コストが高くなります。プロセス数が小さいと、問題が大きすぎてメモリに収まらない可能性があります。こうした理由から、各ベンチマークを異なる解像度と異なる自由度総数で表現して、そのスケール性をプロセス数の最小数の場合に対して相対的に測定します。
そこでテストケースとして以下のデータを選びました：

1. 後方ステップ流れ：このベンチマークでは、2つの異なる解像度を用います。低解像度ベンチマークは250Kノード、高解像度ベンチマークは4Mノードのアダプティブメッシュを用います。アダプティブおよび非アダプティブ版が利用可能です。低解像度ベンチマークはImperial College（IC）とFujitsu Laboratories of Europe Ltd（FLE）の現地リソースで実行されました。高解像度ベンチマークは、RCUKの国立ハイエンドコンピューティングリソースHECToR（Cray XT5）で実行されました。
2. Open ocean deep convection $OODC$ ：360Kノードのアダプティブメッシュを用います。アダプティブおよび非アダプティブ版が利用可能です。
3. 回転円筒：65Kノードのアダプティブメッシュを用います。アダプティブおよび非アダプティブ版が利用可能です。
4. 風成傾圧循環：214Kノードメッシュを64プロセッサー上で動作するように分割します。新しいP1DG-P2要素ペアを用いるため、実行中の実際の自由度はかなり大きくなります。これは約5M速度ノードが存在することを意味します。 多数のコアのスケーラビリティ解析には、200M速度自由度の10M頂点メッシュを用いて大規模ベンチマークを行います。最大で4096プロセスの入力ファイルが生成されます。

ソルバーの比較

ここでは風成傾圧循環データによりソルバーのベンチマークを行います。海洋学における還流とは大規模に回転する海流を指し、これには大規模な風の動きが含まれます。この流れは、コリオリ力と自由表面勾配の水平面のバランスに支配されます。傾圧循環においては、通常は海洋領域の密度の層が考慮されます。これをモデル化するためにFluidity-ICOMは、3D非静水ブシネスク方程式を用います。ベンチマークにおいては如何なる乱流モデルも用いません。
ブシネスク方程式は、圧力要素ペアのP1DGF-P2流速を用いた有限要素積分を通して離散化されます。ここで用いたメッシュは、10M頂点を含み、DG法を用いるため流速自由度は200Mに達します。Fluidity-ICOMは、数値な離散化によるスパース線形系を解くためにPETScライブラリを用います。

還流テストケースでは圧力行列は非常に多くの条件数が存在し、通常のプレコンディショナとソルバーでは手に余る困難さが生じます。ここで、大規模系（その良好なスケール性から）や、大きなアスペクト比を持つ特異な系、または様々な長さスケールを含むより一般的な問題に対するベストな選択としてAMGプレコンディショナを用います。特に大きなアスペクト比を持つ海洋問題を標的にした代数的マルチグリッド法が、Fluidity-ICOMモデルで開発されています$Kramer et al., 2010$。このマルチグリッド法は共役勾配法のプレコンディショナとして用いられます。

PETScはまた、HYPREのBoomerAMG（HYPRE:http://www.llnl.gov/CASC/linear solvers/)のような他のAMG法へのインターフェイスも持ちます。圧力ソルバーに関しては、Fluidity-ICOMのMG法の方がより良好な性能を示します。

ブロック化による行列アセンブリコードの最適化

Fluidity-ICOMの行列アセンブリに対して最適化を行い、アセンブリとベクトル場およびDG場ソルバーに対してblock-CSR格納形式を組み込みました。

メッシュノードの再順序付け

PETScへのインターフェイスで、非零要素が狭いバンド（あるいはクラスタ）になるようにスパース行列の再順序化を行いました。

プロファイリングとスケール性の解析

HECToR上で還流ベンチマークテストケースに対して、CrayPATとVampirを用いたプロファイル採取を行いました。

連続の式による非圧縮性の制約と組み合わせた運動量保存式に焦点を当てます。この計算は、離散化された運動量保存式とポアソン方程式を表現する線形系のアセンブリ、およびこれらの求解から成ります。よって、圧力行列アセンブリと、圧力方程式の線形解法、運動量式アセンブリ、および運動量式の線形解法の4つの部品から構成されます。

1024コアの場合、圧力と流速行列アセンブリは全体の実行時間の30％以上を占め、圧力ソルバーは53.9％に達します。行列アセンブリフェーズが高コストな理由は、深くネストしたループと、その最深ループがサイズに対して４次で増加すること、非構造格子のための関節指標アドレス、およびキャッシュ再利用性によるものです。

1024コアに比較して2048，4096コアでの性能は良好なスケール性を示しています。特に不連続ガラーキン法を用いた流速解法は極めて優れています。圧力アセンブリはこれに比べて効率は下がります。ただしこの処理は全シミュレーションで1回のみ実行されるため、多くの時間ステップを実行する通常のシミュレーションではその影響は僅かです。

·通信オーバーヘッドと負荷バランス

MPI_SYNCは、各集団通信を行うサブルーチンのトレースラッパーで使用され、サブルーチンに入る前にバリア呼び出しで待機する時間が測定されます。 したがって、MPI_SYNCの統計は負荷分散の良い指標となります。
1024から4096までのコア数に対して、MPIで費やされた時間の割合は28.7％から33.1％に増加し、USER関数は45.5％から24.9％に減少し、MPI_SYNCは25.7％から42.0％に増加します。

·高負荷状況

1024コアで最も負荷が高いのは線形ソルバーKSPSolverで、これは4096コアでは3.5倍高速となります。ですが、これは運動量ソルバーの外側ある関数です。
最高負荷のMPI関数はMPI_Allreduceです。この集団通信はスケール性が良くないと予想されましたが、XT4では比較的良好な性能を示します。この関数はPETScライブラリのPetscMaxSumから呼ばれています。MPI_Waitanyは負荷バランスの品質指標で、1024?4096コアの実行でこの量が大幅に増加しないことを考慮すると、コア数の増加において負荷バランスは著しく悪化しているようには見えません。

メッシュアダプティビティのスケール性

アダプティブメッシュ法は、計算領域内のユーザが指定した誤差を達成するメッシュを生成することを目的としています。 ここにはメッシュの辺の長さや、問題のサイズ、ジオメトリ、メッシュのグラデーションに関するユーザ指定の制約が含まれます。 このエラー制御手法は、シミュレーションの計算効率を最適化するものと見ることもできます。 特定のコンピュータ上で静的メッシュを用いて実用的な計算が出来なかった問題を解決することが可能です。

アダプティブメッシュ法を様々な問題に対して並列に処理することが望まれています。メッシュの最適化手法の詳細は文献$Pain et al., 2001$でも見られますが、並列化の大きな課題には、サブ領域間に跨るメッシュの更新や、異なるメッシュ解像度に対するソルバーの負荷バランス、プロセス間のメッシュ移動の最小化という課題があります。

一般にアダプテーションではローカルメッシュ構造の変化を伴います。あるメッシュが更新対象となる場合、そのメッシュに関連する情報を持つ全ての領域が更新されます。階層的メッシュ生成$hierarchical mesh refinement例えばJimack (1998$)とは異なり、このメッシュ最適化方法の性質から、最終的なメッシュのローカルな接続性の異なるプロセスに対する独立した予測は排除されます。

プロセス当たりの相対的な計算負荷の推定は、プロセスに割り当てられたパーティションのメッシュ頂点数から得られますが、これはアセンブルとソルバーの対象となる線形系のサイズに相当します。 並列シミュレーションの冒頭で、グラフ分割法により、プロセス毎のノード数バランスをとる領域分割が行われます。これは、プロセス全体に渡って均一な作業分布を提供するものです。しかしながら、エッジ消去（edge-collapsing $ノード削除$）やエッジ分割（edge-splitting $ノード挿入$）によるメッシュ最適化操作はこのバランスを崩します。こうして引き起こされた負荷分散は数％から数桁のオーダーの差を生じます。このため、グラフ分割法とデータ移動（例えば、新たなグラフ分割を満たすようにプロセスに渡ってメッシュを再分配する）は並列メッシュ最適化における鍵となります。この操作は、メッシュアダプティビティ処理に更なる計算及び通信の負荷を追加します。

このアルゴリズムでは、最初にパーティション上のメッシュは変更されないという制約のもとに、サブ領域は独立にシリアル実行で最適化されます。この最適化が終了した場面では一般に、メッシュの局所的な粗密性が存在するため負荷分散が生じるはずです。よってここで領域の再分割化が必要になります。最適化後の局所的な要素密度の予測にはエラーメトリックフィールド（error metric field）を用いることが出来ます$Pain et al., 2001$。これは、グラフ分割法へ適切なノード重みを定義することでバランスしたパーティション計算を行います。メッシュアダプティビティ適用後は、最適化が適用できないパーティションを除き、これはいずれの場所でも1に近い値になるべきです。

メッシュアダプティビティコードの最適化

以下の手法で最適化します：

グラフ分割手法：
並列処理では、グラフ分割で生成された要素をロックして、各領域のメッシュをアダプテーションします。よってプロセス間の通信は発生しません。

k-wayマルチレベルベース手法：
ここでは、公開されているグラフ分割プログラムParMETIS （ParMETIS: http://glaros.dtc.umn.edu/gkhome/views/metis/parmetis/）を使用しました。 これは、グラフのエッジと頂点の両方にウェイトを適用可能な高速k-wayマルチレベルベース手法を使用します。 ここでの議論に関連する3つの主要な並列再分割方法は、scratch、scratch-remappingおよびdiffusion手法です（Schloegel et al., 1997）。

エッジ重みとノード重みの最適化：
グラフ分割の適切な摂動は、再分割時にグラフにエッジ重みを適用することによって達成されます。以前に（グラフ分割のために）ロックされた領域に高いエッジ重みを適用するだけで、以前にロックされた領域のすべてが最適化を必要としないため、プロセス間の不要なメッシュ移行を促進することができます。

データ移動：
再分割後に必要な部分をプロセス間で移動させます。

既存のパーティションと新しいパーティションの照合：
通信に必要な最小のデータセットを決定します。

通信の最適化：
移動すべき情報をシリアライズ化して通信します。

パーティション再構築：
ハロ情報と共に新たなパーティションを構築します。

メッシュ最適化アルゴリズムの評価

ここでは強スケール性に着目して、実計算における並列メッシュ最適化手法の性能評価を行います。

差動的加熱回転円筒ベンチマークケースを用いて、次の2つの異なるサイズを検討します。これは、1, 2, 4, 8, 16, 32, 64, 128プロセッサを用いた0.1Mノードメッシュ、および64, 128, 256, 512, 1024プロセッサを用いた2.0Mノードメッシュのケースです。この測定では、運動量と温度計算のソルバーとプレコンディショナとして、PETScのGMRESとEISENSTATを用いました。圧力ソルバーとしては、共役勾配法とHYPRE's BoomerAMGプレコンディショナを用いました。

これらの実験では、メッシュアダプティビティ処理は全実行時間の約5％を占めています。この時間はシミュレーションの問題そのものの特性と、アダプティビティの適用頻度に依存します。全ての実験においてプロセス数の増加と共に高速化が観測されました。

プロセス数の増加に対する全体的な並列効率の低下は、問題サイズを固定した場合の、プロセス数が増えるにつれて計算/通信比が減少することと対応しています。これは、パーティションが所有するノード平均数とノード総数との相対的な比として見ることが出来ます。たとえば、128プロセスで100kノードの問題を実行すると、パーティションあたり約800個のパーティション所有ノードと約700個のハローノードが作成されます。 1024プロセスで2Mノードの問題を実行すると、パーティションあたり約2000のパーティション所有ノードと約1400のハローノードが発生します。調査により、効率の低下は、シリアル実行のメッシュアダプティビティ操作の不均衡さでなく、データ移動/再バランス処理のオーバーヘッドによるものであることが明らかになりました。このベンチマークは、実行されるメッシュアダプティビティ操作がドメイン全体でほぼ等配分であったため、最大限のデータ移行量が必要な点で限定的なケースと見なすことができます。

Vampirによりさらに詳細なプロファイリングを行い、ParMETISで複数の冗長なゼロサイズメッセージ通信（MPI_IsendとMPI_Ireceiveを使用しています）が明らかになりました。アダプティビティにおけるデータ移動コストは非常に低いため、現時点ではさらなる調査を行いませんでした。

並列I/O

出力ファイルはVTK$http://www.vtk.org/$で管理されています。VTKは非構造格子に対して.vtuファイルフォーマットを適用します。N回のダンプ回数に対して、simulationname_N.vtuという出力結果を含むファイルが生成されます。このファイルには、ダンプ時間直前の実行のスナップショットまでが含まれます。例えば、時間ステップが3秒に設定され、10秒毎にダンプ指定がされている場合、最初のダンプは12秒後です。

並列シミュレーション時には、データは.vtuおよび.pvtuファイルの両方に保存されます。.vtuファイルには各パーティションの出力データが含まれ、ファイル名はsimulationName_P_N.vtuです（Pはプロセッサ番号、Nはダンプ番号です）。 .pvtuファイルは（ダンプ番号毎の）全メッシュに関する出力データを含みます。

.vtuファイル出力はインターリーブI/Oを行います。Nコア（Nプロセス）を用いる場合、各プロセスは複数のローカルファイルを出力し、全体で5N（3Dメッシュの場合）ファイル生成されます。この場合、LUSTREのような並列ファイルシステムによるシンプルな並列I/O実装が期待できます。しかしながら、大規模プロセス数（ファイル数）ではメタデータ操作、さらにOSSおよびOSTの競合が全体性能を阻害します。

インターリーブI/Oは、アプリケーションプロセスのサブセットを使用してI/Oを実行します。 これにより、I/O段階でメタデータサーバーに到達する要求の数が制限され、パフォーマンスが向上します。 ライターの数とランクオーダーの両方を慎重にチューニングすると、パフォーマンスが向上します。

還流ケース（10Mメッシュ頂点）で解析した結論として、4096プロセスまでは、デフォルトのI/OがインターリーブI/Oに勝ることが判りました。

謝辞

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

文献

[1]	Cotter, C., Ham, D., Pain, C., 2009. A mixed discontinuous/continuous finite element pair for shallow-water ocean modelling. Ocean Modelling 26, 86-90.
[2]	Ford, R., Pain, C., Piggott, M., Goddard, A., de Oliveira, C., Umpleby, A., 2004. A Nonhydrostatic Finite-Element Model for Three-Dimensional Stratified Oceanic Flows. Part I: Model Formulation. Monthly Weather Review 132 $12$, 2816-2831.
[3]	Gorman, G., Pain, C., Piggott, M., Umpleby, A., Farrell, P., Maddison, J., 2009. Interleaved parallel tetrahedral mesh optimisation and dynamic load-balancing. Adaptive Modelling and Simulation 2009, 101-104.
[4]	Gorman, G., Piggott, M., Wells, M., et al, 2008. A systematic approach to unstructured mesh generation for ocean modelling using GMT and Terreno. Computers & Geosciences 34, 1721-1731.
[5]	Gorman, G. J., Piggott, M. D., Pain, C. C., de Oliveira, C. R. E., Umpleby, A. P., Goddard, A. J. H., 2006. Optimisation based bathymetry approximation through constrained unstructured mesh adaptivity. Ocean Modelling 12$3-4$, 436-452.
[6]	Gresho, P., Sani, R., 1998. Incompressible Flow and the Finite Element Method.
[7]	Hendrickson, B., Leland, R., 1995. A multilevel algorithm for partitioning graphs. In: Proc. Supercomputing '95. Formerly, Technical Report SAND93-1301 $1993$.
[8]	Jimack, P. K., 1998. Techniques for parallel adaptivity. In: Topping, B. H. V. $Ed.$, Parallel and Distributed Processing for Computational Mechanics II. Saxe-Coburg Publications.
[9]	Karypis, G., Kumar, V., 1998. A fast and high quality multilevel scheme for partitioning irregular graphs. SIAM Journal on Scientific Computing 20 $1$, 359-392.
[10]	Kernighgan, B. W., Lin, S., 1970. An efficient procedure for partitioning graphs. Bell Systems Technical Journal 49.
[11]	Kramer, S., Cotter, C., Pain, C., 2010. Solving the Poisson equation on small aspect ratio domains using unstructured meshes. Ocean Modelling, submitted, arXiv:0912.2194.
[12]	Kramer, S., Pain, C., 2010. Modelling the free surface using fully unstructured meshes. In preparation.
[13]	Ozturan, C., 1995. Distributed environment and load balancing for adaptive unstructured meshes. Ph.D. thesis, Rensselaer Polytechnic Institute, Troy, New York.
[14]	Pain, C., Piggott, M., et al., 2005. Three-dimensional unstructured mesh ocean modelling. Ocean Modelling 10, 5-33.
[15]	Pain, C., Umpleby, A., de Oliveira, C., Goddard, A., 2001. Tetrahedral mesh optimisation and adaptivity for steady-state and transient finite element calculations. Computational Methods in Applied Mechanics and Engineering 190, 3771-3796.
[16]	Schloegel, K., Karypis, G., Kumar, V., 1997. Multilevel diffusion schemes for repartitioning of adaptive meshes. Journal of Parallel and Distributed Computing 47 $2$, 109-124.
[17]	Selwood, P. M., Berzins, M., 1999. Parallel unstructured tetrahedral mesh adaptation: algorithms, implementation and scalability. Concurrency: Practice and Experience 11 $14$, 863-884.
[18]	Walshaw, C., Cross, M., Everett, M., 1997. Parallel dynamic graph partitioning for adaptive unstructured meshes. Journal of Parallel and distributed Computing 47, 102-108.