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

チューニングレポート：大規模並列アーキテクチャにおける雲/エアロゾル・シミュレーションの並列化

Jon Gibson, Numerical Algorithms Group Ltd
28 February, 2013


[2016年6月掲載]

概要

本プロジェクトはHECToR：英国国立学術スーパーコンピューティングサービスにおける1年のサポートプロジェクトとして遂行されました。本レポートは最初にイギリス気象庁(Met office(MO))Large-eddyモデル(LEM)内への並列繰返し法ソルバーの実装作業を記述します。この新しいアプローチはFFTベースのソルバーの代替案として提案され、側面の周期境界条件を除いて新たな科学計算を可能にするものです。次に、エアロゾルからの雲の生成と成長を研究するためのプロセススケールモデルである、コードACPIM(Aerosol-Cloud-Precipitation Interaction Model)の並列化について記述します。プロジェクト責任者はリーズ大学のAlan Gadian博士とマンチェスター大学のPaul Connolly博士です。

目次



1.　アプリケーションの概要：LEM

Large Eddy Model(LEM)コードはイギリス気象庁で開発され、エアロゾルと雲の相互作用解析と言った雲スケールおよびミクロスケールの大気プロセスを高解像度数値モデルを用いてシミュレートします。乱流、微視的物理プロセスおよび放射をパラメータ化した修正Navier-Stokes方程式が実装されています。このソフトウェアはイギリス気象庁外の科学コミュニティでも幅広く利用されており、文献[2,3,4]に文書化されています。

1.1　基礎方程式系

文献[3]において詳細に記述されていますが、LEMは以下の基礎方程式系を解きます。ここでテンソル記法を用います。これら方程式は中心差分スキームを用いて時間積分されます。

\begin{eqnarray} \frac{Du_i}{Dt} = -\frac{\partial}{\partial x_i} \left( \frac{p'}{\rho_s} \right) +\delta_{ij}B'+\frac{1}{\rho_s}\frac{\partial \tau_{ij}}{\partial x_j} -2\epsilon_{ijk}\Omega_ju_k \tag{1} \\ \frac{\partial}{\partial x_i}(\rho_s u_i) = 0 \tag{2} \\ \frac{D\theta}{Dt} = \frac{1}{\rho_s}\frac{\partial h_i^{\theta}}{x_i} +\left( \frac{\partial \theta}{\partial t} \right)_{mphys} +\left( \frac{\partial \theta}{\partial t} \right)_{rad} \tag{3} \\ \frac{Dq_n}{Dt} = \frac{1}{\rho_s}\frac{\partial h_i^{q_n}}{x_i} -\left( \frac{\partial q_n}{\partial t} \right)_{mphys} \tag{4} \end{eqnarray}

ここで、

\begin{eqnarray} \boldsymbol{u} &:&　流速ベクトル \\ \theta &:&　温位 \\ q_n &:&　他の全てのスカラー変数の代表 \\ p &:&　圧力 \\ \rho &:&　密度 \\ B’ &:&　浮力 \\ \tau &:&　サブグリッドの応力 \\ h^{\theta} &:&　サブグリッドの\thetaのスカラーフラックス \\ h^{q_n} &:&　サブグリッドのq_nのスカラーフラックス \\ \Omega &:&　地球の角速度(f平面近似) \\ \epsilon_{ijk} &:&　交代擬テンソル \\ \left( \frac{\partial \theta}{\partial t} \right)_{mphys} &:&　微視的物理による\thetaのソース項 \\ \left( \frac{\partial \theta}{\partial t} \right)_{rad} &:&　放射による\thetaのソース項 \\ \left( \frac{\partial q_n}{\partial t} \right)_{mphys} &:&　微視的物理によるq_nのソース項 \end{eqnarray}

および以下を用いた。

\begin{eqnarray} \frac{D}{Dt} \equiv \frac{\partial}{\partial t}+u_i \frac{\partial}{\partial x_i} \tag{5} \end{eqnarray}

1.2　圧力計算

式(1)の圧力項の計算は下記のポアソン様楕円型方程式の求解が必要です。

\begin{eqnarray} \frac{\partial}{\partial x_i} \left[ \rho_s \frac{\partial}{\partial x_i}(p’/\rho_s) \right] = \frac{\partial}{\partial x_i}(\rho_s s_i) \tag{6} \\ ここで、\\ s_i=\delta_{ij}B’-u_j \frac{\partial u_i}{\partial u_j} +\frac{1}{\rho_s} \frac{\partial \tau_{ij}}{\partial x_j} -2\epsilon_{ijk} \Omega_j u_k \tag{7} \end{eqnarray}

これは水平面でフーリエ変換をしてx,yの偏微分を消去し、その結果を三重対角ソルバーで解き、最後に逆フーリエ変換をしてp’/ρsを得ます。この方法がプロジェクト前にコードに実装されていた唯一の方法です。

1.3　LEMソフトウェア

このソフトウェアは文献[4]に詳述されていますが、ここでコード特性に適切な文脈において記載することとします。このコードは主にFortran77で記述され、一部がFortran90に従っています。全メモリーは静的に確保されています。メッセージ通信は、MPIライブラリの小さなサブセットのラッパー(MPLと呼ばれます)のラッパーであるGCOMライブラリ[5,6]を用いて記述されています。このアイデアは、低レベルのMPI[7]、PVM[8]、SHMEM[9]を同じGCOM呼び出しで可能にしたものです。しかしながら現在MPIが広く使われるようになり、GCOMはMPI機能の小さな部分のみを用いているため、効果的なプログラミングの障害になっています(MPI2のみならずMPI1の多くの機能も持ち合わせていません)。

バージョンコントロールは古い手法を取っています。用いているのはnupdate source code manager[10]です。このシステムで扱うソースコードはdeck(,dkファイル)とcomdeck(.cdkファイル)から成り、nupdateコマンド実行によりFortranソースコードを生成させます。変数とcommon宣言を含むcomdeckは、その処理中にソースコードへインライン展開されます。ソースコードへの如何なる修正も、他のファイルとともに処理されるupdateファイル書き出しが必要です。

実行jobはlemsubスクリプトを用いてユーザーによりsubmitされます。これはnupdateコマンド、デフォルト情報ファイル、マシン情報ファイル、マシン依存実行スクリプト、makefile、updateを通して、実行ディレクトリへのファイル転送、コンパイル、そしてバッチシステムへのjob投入と続きます。この初期jobをset-up-jobと呼ぶこととします。引き続きjob投入を行うように準備される場合はchain-jobと呼ばれ、特定のポイントまでシミュレーションが到達するまで繰り返されます。

1.4　並列領域分割

問題領域はIIP×JJP×KKPのグリッドポイントで準備されます。LEMはy-z平面の2Dスライスをx方向へ連続して動かして場を積分していきます。モデル方程式の積分には高々2平面の情報しか要求しないので、領域はオーバーラップするharo領域を持つサブ領域へ分割可能です。図1は文献[4]から取ったもので、問題領域がx方向にどのように分割されるかを図示しています。各PEはIIP/NPES枚のスライスを持ち、隣接サブ領域には2つのスライスに対するハロ領域を持っています。


図1:LEMの領域分割。図aはシングルプロセッサ実行を示す。ハロ領域スライスは周期境界条件を満たす(スライス1はスライスIIP-1のコピー、スライス0はスライスIIPのコピー、等)。図bは4プロセッサでの実行で、ハロ領域は隣接プロセッサ間の連続性を保つために用いられる(PE0上のIIPEP+1とIIPEP+2は、PE1上のスライス1,2のコピーである)。

2　アプリケーション概要：ACPIM

エアロゾル-雲-降水相互作用モデル(ACPIM)はモジュラーフレームワークを持ち、マンチェスター大学、およびドイツの気象・気候研究所における世界最先端のAIDA雲チャンバ施設で、開発および検証されています。ACPIMは、エアロゾルからの雲の生成と成長を研究する詳細なプロセスモデルとして主に用いられてきました。このコードはC言語で記述されたシリアルコードです。

3　プロジェクトの概要

以下の目標が最初に提案されました。

目標1：LEMの圧力ソルバーを、HECToRに適した並列FFTライブラリーを用いて高速化する。
ポアソン方程式を解くためには最初にFFTを用い、ポアソン方程式に現れる2つの水平方向の微分を消去して、水平グリッドの各ポイントでの2階のODEへ変換します。このODEを、最初に3つの鉛直グリッドポイントで離散化して、次に三重対角ソルバーを用いて逆行列により解きます。最後に逆フーリエ変換により圧力の周波数モードを大気圧力へ変換します。現状のLEMはシリアルFFTソルバーを用いていますが、本プロジェクトではHECToRへ最適化するためにFFTWのような並列バージョンを実装します。この手法は高速ですが、周期解しか解けないため場合によっては使いずらいものです。繰返し法ソルバーは周期境界条件は不要です。

目標2：LEMへHECToRに適した並列繰返し法ソルバーを実装する。
別のポアソン方程式解法として繰返し法があります。GKRES法,GMRES法が効率的であることが知られており、周期境界条件も不要です。米国ではマルチプロセッサ・アーキテクチャの利用において、ポアソン方程式解法に繰返し法ソルバーが必要であることが認められています。大気流の非弾性方程式に対するGCR(k)マルチグリッドソルバーは、Thomas & Smolarkiewicz (1997, Copper Mountain Conference on Multigrid Methods, http://www.mgnet.org/mgnet-ccmm97.html)により最初に提案されました。これは非自己共役作用素を用いた方程式の非振動性の時間方向前進差分(forward-in-time (NFT))による離散化に基づいており(Smolarkiewicz & Margolin 1997, Atmos Ocean, 35,127-157)、スペクトル法プレコンディショナが効率を改善することが示されています(Thomas SJ, et al Mon. Wea. Rev 131(10): 2464)。本プロジェクトでの提案はGMRES法(Allen 2007, Atmos Sci. Letters, DOI , 10.1002/asl.145)です。GMRES法の利点は、すでに実装されベクトル演算のみですが稼働しているUK MO Blasius modelと同等であることです。ここでもプレコンディショナが効率向上のため実装されています。完成後にイギリス気象庁の研究者による発表が予定されています。この手法はイギリス気象庁で利用可能になる予定で、採用されればこのLEMモデルの他のユーザにも利用可能になるでしょう。

FFTおよび繰返し法ソルバー実装には9人月のサポートが必要とされました。このほとんどは繰返し法の実装で占められるでしょう。

目標3： HECToR向けにACPIMモデルを並列化する。
ここでの主要な課題は、良好なスケーラビリティを示すようにMPIを実装することです。

4　並列繰返し法ソルバー

本プロジェクトの最も大きな作業で、FFTベースの手法に代わり、圧力解法のために並列繰返し法ソルバーを実装します。この手法の主要な利点は周期境界条件が不要なことで、LEMの能力を拡張するものです。繰返し法ソルバーはポアソン方程式形式の式を解きます。

\begin{eqnarray} \nabla^2 \boldsymbol{x}=\boldsymbol{b} \tag{8} \end{eqnarray}

多くのやり方の中で3種の方法を実装しました：BiCGStabとGMRESおよび、GCRです。その比較から明らかに最も効率が良いのはBiCGStab(BiConjugate Gradient Stabilized)法でした。これはイギリス気象庁の示唆[11]とも一致します。

4.1　BiCGStabアルゴリズム

ここで行列Aを用いて、Ax=∇2x=bとして書くこととします。アルゴリズムは以下の通りです。

  1. 1. 終結交差TOLと最大繰返し数ITMAXをセットする。
  2. 2. 誤差に対するスケール因子$K_{err}=max(\sqrt{\boldsymbol{b} \cdot \boldsymbol {b}},0.01)$を求める。
  3. 3. 初期値$\boldsymbol{x}_0$に対する初期残差$\boldsymbol{r}_0=\boldsymbol{b}-\boldsymbol{Ax}_0$を求める。
  4. 4. 初期化：$\bar{\boldsymbol{r}}=\boldsymbol{r}_0, \boldsymbol{p}_0=\boldsymbol{\nu}_0=0,\rho_{-1}=\alpha_0=\omega_0=1$
  5. 5. DO i=1,ITMAX
  6.   a. $\rho_{i-1}=\bar{\boldsymbol{r}}^T \cdot \boldsymbol{r}_{i-1}$
  7.   b. $\beta_{i-1}=\left( \frac{\rho_{i-1}}{\rho_{i-2}} \right)\left( \frac{\alpha_{i-1}}{\omega_{i-1}} \right)$
  8.   c. $\boldsymbol{p}_i=\boldsymbol{r}_{i-1}+\beta_{i-1} [\boldsymbol{p}_{i-1}-\omega_{i-1} \boldsymbol{\nu}_{i-1}]$
  9.   d. $\boldsymbol{M \hat{p}}=\boldsymbol{p}_i$を解く
  10.   e. $\boldsymbol{\nu}_i=\boldsymbol{A \hat{p}}$
  11.   f. $\alpha_i=\frac{\rho_{i-1}}{\boldsymbol{\hat{r}}^T \cdot \boldsymbol{\nu}_i}$
  12.   g. $\boldsymbol{s}=\boldsymbol{r}_{i-1}-\alpha_i \boldsymbol{\nu}_i$
  13.   h. $\boldsymbol{M \hat{s}}=\boldsymbol{s}$を解く
  14.   i. $\boldsymbol{t}=\boldsymbol{A \hat{s}}$
  15.   j. $\omega_i=\frac{\boldsymbol{t}^T \cdot \boldsymbol{\hat{s}}}{\boldsymbol{t}^T \cdot \boldsymbol{t}}$
  16.   k. $\boldsymbol{x}_i=\boldsymbol{x}_{i-1}+\alpha_i \boldsymbol{\hat{p}}+\omega_i \boldsymbol{\hat{s}}$
  17.   l. $\boldsymbol{r}_i=\boldsymbol{\hat{s}} -\omega_i \boldsymbol{\hat{t}}$
  18.   m. 残差の収束判定：IF( $\frac{\sqrt{\boldsymbol{r}_i \cdot \boldsymbol{r}_i }}{K_{err}} \lt TOL$)  THEN   exit   loop
  19. 6. END   DO

ステップ3,5e,5iでの$\boldsymbol{Af}$の計算では、7点ステンシルを用いた以下の有限差分式を仮定しています：

\begin{eqnarray} \boldsymbol{Af}=\nabla^2\boldsymbol{f}= \frac{1}{\Delta x^2}\left[ f(x+\Delta x,y,z)+f(x-\Delta x,y,z) \right] \\ +\frac{1}{\Delta y^2}\left[ f(x,y+\Delta y,z)+f(x,y-\Delta y,z) \right] \\ +\frac{1}{\Delta z^2}\left[ f(x,y,z+\Delta z)+f(x,y,z-\Delta z) \right] \\ -2\left( \frac{1}{\Delta x^2}+\frac{1}{\Delta y^2}+\frac{1}{\Delta z^2} \right)f(x,y,z) \tag{9} \end{eqnarray}

ここで、$\Delta x,\Delta y,\Delta z$はそれぞれz,y,z方向の近傍グリッド点への距離を表します。LEMではグリッド間隔はx,y,z方向に一定ですが、z方向に動くのに従い変化します。

ステップ5dと5hはプリコンディショナを用います。ヤコビ法では以下に従います[11]：

\begin{eqnarray} x_i^{(k+1)}=\frac{1}{a_{ii}} \left( b_i- \displaystyle \sum_{ j \neq i }^{}a_{ij}x_j^{(k)} \right) \tag{10} \end{eqnarray}

このコードでは、プレコンディショナの繰返し回数はパラメータPREITSで指定されます。

4.2　並列BiCGStab

並列BiCGStabにおいて、全てのLEMユーザはx方向で同じ領域分割を用いることを仮定します。これが意味するのは、式9の$\boldsymbol{Af}$の計算が隣接プロセスとハロ交換を要求することです。ステップ2,5a,5f,5j,5mには内積計算が含まれますが、これは全プロセスに渡る全総和計算を要求します。これ以外の演算は全てローカルな計算です。

代替手法も可能ですが、当面はx,y方向には周期境界条件が想定されています。z方向の圧力境界値は0にセットされます。

4.3　実装

LEMバージョン2.4へ繰返し法ソルバーが追加されました。以下のファイルがソースコードへ追加されました。その全てのdeckファイルは、.dkを除いた名前の単独のプロシージャを含みます。これらは、関数INNER_PROD以外は全てサブルーチンです。

  1. POISSON_ITERATIVE.dk：　繰返し法ソルバーを呼び出して圧力ポアソン方程式を解く。これはFFTでポアソン方程式を解くPOISSON.dkと同等機能である。
  2. SET_MAT_A.dk：　Ax=∇2xにおける行列Aをセットする。
  3. BICGSTAB.dk：　BiCGStab法による繰返し法ソルバー。実験の結果から、パラメータITMAX,TOL,PREITSをそれぞれ999,10-8,1にセットした。ユーザがこの値を変えることは想定してないが、必要となればこのファイルを編集して変更できる。
  4. PRECOND.dk：　ヤコビ法プレコンディショナ。PREITS回実行する。
  5. CALC_Ax.dk：　Ax=∇2xを計算する。
  6. INNER_PROD.dk：　全プロセスに分散している2配列の内積を計算する。
  7. SOLVER.cdk：　繰返し法ソルバーに関係する変数を格納するcomdeckファイル。

次のソースコードファイルは繰返し法ソルバー導入前からあるものですが以下のように更新されています：

  1. NNSTEPS.dk：　プリプロセッサーで解釈される#ifdef　#else　#endif文が挿入された。マクロ：ITERATIVE_SOLVERが定義された場合はサブルーチンPOISSON_ITERATIVEが呼び出され、そうでなければ標準バージョンとしてPOISSONが呼び出される。

以下のファイルもソースコードディレクトリに存在するものですが、追加のファイルを追加して更新されました：

  1. deck.list：　nupdateで処理されるdeckファイルリスト。
  2. comdeck.list：　nupdateで処理されるcomdeckファイルリスト。

図2:並列繰返しソルバーで用いるルーチンのcallグラフ。

4.4　コード検証

繰返し法ソルバーの正しさは、テストコードを用いて以下のように確認されました。

  1. 1. 既知の圧力値φを作成する
  2. 2. f=∇2φにより入力ソース項を計算する。
  3. 3. fを入力として繰返し法ソルバーを呼び出す。
  4. 4. この結果が元の圧力値φと一致するかをチェックする。

4.5　コード最適化

Cray PATプロファイルツールとコンパイラ出力結果を用いてコードの性能改善、主にキャッシュ効率の改善を行いました。繰返し法ソルバーの通信は十分に効率的でしたが、それ以外のコード部分は効率が悪く、修正が必要でした。プロファイリング作業の中で、ノンブロッキング通信がハロ交換の性能を改善することが明らかになりました。GCOMにはこの基本的なMPIの機能が含まれないため、必要に応じて中間的なMPLライブラリを直接呼び出すように修正しました。

4.6　コード性能

図3にHECToRフェーズ1上での繰返し法ソルバーによるLEMのスケール性能を示します。


図3:　 MPIプロセス数に対する、HECToRフェーズ1上での繰返し法ソルバーによるLEMの実行時間[秒]のグラフ。線形にスケールしている。

5.　並列ACPIM

コードは2重のメインループを持ちます。外側は時間ステップ、内側は1D列を移動する空気塊(air parcel)のループです。このparcelセットはMPIプロセスで分割されます。Parcel間には通信は発生しませんが、そのデータは時間ステップ毎にルートプロセスへ集められ、そのデータをファイルへ出力し、遺留と沈降計算を行う必要があります。並列効率はシミュレーションで用いるparcel数で制約され、通常は20から30程度です。

6.　結論

MO LEMコードへ成功裏に繰返し法ソルバーが追加されました。これは以前のFFTを用いる場合に必要な周期境界条件を不要にしました。これにより新たな研究領域において利用が可能となりました。

ACPIMコードも並列化され、同様に新たな研究領域への道を開きました。

3つ目のプロジェクト目標はLEMコード内のFFTをFFTWライブラリ呼び出しに変更することでしたが、時間制約とGNUライセンスから、イギリス気象庁のコードに採用することは適切でないと判断されました。

謝辞

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

文献

[1] http://www.hector.ac.uk/cse/
[2] M. E. B. Gray and J. Petch, “Version 2.3 Of The Met Office Large Eddy Model: Part I. User Documentation”
[3] M. E. B. Gray et al, “Version 2.3 Of The Met Office Large Eddy Model: Part II. Scientific Documentation”
[4] M. E. B. Gray et al, “Version 2.3 Of The Met Office Large Eddy Model: Part III. Software Documentation”
[5] Jørn Amundsen and Roar Skålin, “GC User’s Guide Release 1.1”, SINTEF Applied Mathematics, October 16, 1996.
[6] Jørn Amundsen and Roar Skålin, “GCG User’s Guide Release 1.1”, SINTEF Applied Mathematics, October 16, 1996.
[7] http://www.mpi-forum.org/docs/docs.html
[8] http://www.csm.ornl.gov/pvm/
[9] http://staff.psc.edu/oneal/compaq/ShmemMan.pdf
[10] Ramesh Krishna, “Unified Model Documentation Paper X2: The Nupdate Source Code Manager”, Met Office, November 11, 1993.
[11] Meeting with Tom Allen, Met Office.
関連情報
MENU
Privacy Policy  /  Trademarks