＊ここに掲載するのは、サザンプトン大学のR. Johnstone博士によるHECToRレポート「Improved Scaling for Direct Numerical Simulations of Turbulence, Johnstone, R. University of Southampton Southampton, SO17 1BJ, August 2, 2012」を要約したものです。
概要
ここで取り上げるのは2つのプログラムです。一つはSWTで、無限平面内の乱流の直接数値シミュレーション用に設計されました。これは最近になって、非回転ひずみに対する流れのレスポンス[3]や乱流クエットポアズイユ流（移動壁や圧力勾配両方から生起されるチャネル流[9]）のシミュレーションにも用いられました。これらは乱流上の圧力勾配の影響を明らかにする研究に役立てられます。第2のコードはSS3Fで、これまで主に成層流内の渦ダイナミクスのシミュレーションに用いられています[1,2,4]。
SWTはKim et al [10]のベクトル性[15]チャネル流コードで、並列化された現代的なコードです[3]。SWTはSS3Fと同じルンゲクッタ法とFFTルーチンを用います。一方向にはチェビシェフ基底を用いますが、一方向のみ非周期的な平面チャネル流に適しています。
SS3Fはブジネスク近似の非圧縮性ナビエストークス方程式を解き、3次元全ての方向にフーリエ表現を用います。Corral & Jiménez [5]による無限境界条件がこれらの内の一つに実装され、その他は周期的に扱われます。粘性項は解析的に時間進行し、低容量の3次のルンゲクッタスキームが移流項と浮力項に用いられます。
数値計算方法
両方のコードは共に、N-S方程式の非線形移流項を擬スペクトル法で解き、波数空間との間の変換はスライス（1次元データ分割）内で実行します。これは、残りの方向に関する変換には並列転置（all to all通信）が要求されることを意味します。
ナビエストークス方程式は4次の非線形のため、ゼロで無い係数は波数±2kmaxまでです。ここで元の波数範囲は[-kmax:+kmax]です。非線形項を波数空間へ変換する際に、この範囲外の係数を範囲内へ混入するのを避けるようにします。
並列化
このプロジェクトの前は、SWTとSS3Fの両方は1D領域分割を用いていました。領域は2Dスライスに分割されてFFTは2次元で実行されるため、別の一次元のみ異なる2Dスライスへの並列転置操作が必要です。
ここでは3つの連続した1D変換として3Dフーリエ変換を実装しました。つまり個々のプロセッサは一度に1次元のみ扱えば良いことになります。領域を2つの不完全な次元を持つ要素である"ペンシル"へ分割する2D領域分割は、明らかに利用可能なプロセッサ数を増やします。
しかしながらこの手法はコスト的影響があります。変換の度に並列転置操作が必要になります。よって、1D変換は少数プロセッサ数では好ましい手法として残りますが、プロセッサ数は高々、1次元方向のサンプルポイント数に制限されます。またグリッドが細かくなればメモリーや計算負荷はグリッド数の2乗で増加します。一方プロセス数は線形にスケールするのみです。
実際には状況はより好ましくなく、一般に2つの異なる1D分割が必要で、もし各方向のグリッドポイント数が同じでないなら、プロセッサー数は小さい方が選択されます。さらに、その何方かがプロセッサ数で割り切れない場合は負荷分散が生じ、こちらの方が問題になります。
これらのシミュレーションに対する制約から、SWTとSS3Fは現状では数百プロセッサコアに制限され、シミュレーション実行時間的に低い性能しか持ち得ません。より深刻なのはメモリーのスケール性です：全3次元方向の解像度を2倍にするとコア数は2倍にできますが、全メモリーは8倍になります。マルチコアアーキテクチャの流行からコア当たりのメモリー比の減少を期待するだけでは問題は悪化するのみです。
乱流シミュレーションの現実世界への適用には、グリッド解像度が決定的です（反例として[14]を参照してください）。そうでなければ、忠実度が必須であるか、対象とする問題を修正する：例えばレイノルズ数を下げるなどが必要です。
チャネル流の場合SWTは、散逸レートεが壁における摩擦速度uτとチャネル高hを用いて仮定をしています。
この仮定から、シミュレーションに要求されるメモリーはRe9/4に比例することが導かれます。最小の渦を解像するには時間ステップも小さくしなくてはなりません。
同様な仮定をして、最大時間スケール：h/平均流速時間（uτに比例すると仮定）と、最小時間スケールであるコルモゴロフ時間スケールの時間比は、Re1/2に比例することが導かれます。これはKim et al [10]による結果と同じです。
これはレイノルズ数に関する計算負荷のスケール性、つまり2乗に比例することを示しています。
作業内容
このプロジェクトは2011年11月に開始され、HECToRフェーズ2b上で性能分析を行いました。
既存コード性能
既存コードは並列転置をノンブロッキングMPI送受信で実装しています。Crayマシン上ではMPI_ALLTOALLV（2DECOMP&FFTが使用）の方が性能が良いため、置き換えが効果が在ると期待できます。
325コア上での864 x 325 x 324グリッドを用いたSWTコードのプロファイリング結果から、実行時間の内、MPI_WAITに25％、MPI_SYNCに10％掛かっています。これは負荷分散のために生じている可能性が大きく、2D分割を導入すれば改善すると考えられます。
プロファイルのユーザセクションにおいて、殆どがFFTパッケージ(vecfft)であり、その効率化が重要です。最も負荷の高いFFTルーチンがrad3bですが、キャッシュ効率が90.1%と低く、SSE命令によるベクトル化が全くされていません。
作業計画
既存コードのプロファイルからSWTとSS3Fは、FFTの最新化と2D領域分割によりコア当たりの性能と並列スケール性の改善が期待出来ることがわかりました。
既存FFT(vecfft)のFFTW3への置き換えは、シンプルに個別の一次元変換へ置き換えることとしました。チェビシェフ変換のための多重コサイン変換も試しましたが性能改善を示しませんでした。
SS3Fの実行時間はこの置き換えの結果、256 x 1440 x 2880サイズで55％削減されました。
SWTでは3072 x 325 x 1024サイズで53％削減されました。
SWTにFFTWを用いる様に変換すると、チェビシェフ変換ルーチンが必要になります。2DECOMP&FFTのFFTのAPIを用いるには、これをライブラリーに組み込む必要があります。一方領域分割APIを使う場合にはそれは必要ありません。しかしながら組み込んだ方が今後の展開に便利なため組込を行いました。
チェビシェフ多項式計算にはChebyshev-GaussおよびGauss-Lobatto点を実装し、ユーザが選択出来る様にしました。
頂点の再番号付与
良好なメッシュ分割の後には、2つやるべきことが有ります。最初に、新しい分割に応じた要素データの再配分を行います。このため、各プロセスにおいて、他の全てのプロセスに配るべき要素を格納したバッファをパックします。このバッファはMPI_Alltoallvで交換され、その後バッファを新しい要素データ構造へアンパックします。
次に、新しい要素分割に応じた頂点を再度番号付けする必要があります。これはプロセス内ローカルな頂点番号を連続にするためです。プロセス数の増加に伴いインターフェイス要素が重要になります。この時表面/体積比は増加し、インターフェイスの頂点の所属を決めるやり方が重要になります。
ここでシリアル版と並列版の再番号付けアルゴリズムを実装しました。並列版は150百万個の要素ケースでは、プロセス数に依らず20秒で完了ました。
結果
SS3Fは、12000コアまで良好にスケールしました。18000コアでは若干性能が劣化しています。
SWTは、256ノード（8129コア）迄行くとスケール性が劣化しますが、それでもコストは34％少なく、言い換えれば30倍高速と言えます。
謝辞
このプロジェクトは、nAG Ltd.が運営するHECToRの分散計算科学および工学(CSE)サービスの基に実行されました。英国の国立スーパーコンピューティング・サービスである、HECToR：英国リサーチ・カウンシル・ハイエンド計算サービスは、リサーチ・カウンシルを代行するEPSRCが管理しています。そのミッションは英国学術界の科学および工学の研究支援です。HECToRスーパーコンピューターは、UoE HPCx Ltd.およびnAG Ltd.のCSEサポートサービスにより管理運営されています。
文献
| [1] | P.J. Archer, T.G. Thomas, and G.N. Coleman. Direct numerical simulation of vortex ring evolution from the laminar to the early turbulent regime. Journal of Fluid Mechanics, 598:201-226, 2008. |
| [2] | P.J. Archer, T.G. Thomas, and G.N. Coleman. The instability of a vortex ring impinging on a free surface. Journal of Fluid Mechanics, 642:79-94, 2010. |
| [3] | G.N. Coleman, D. Fedorov, P.R. Spalart, and J. Kim. A numerical study of laterally strained wall-bounded turbulence. Journal of Fluid Mechanics, 639:443-478, 2009. |
| [4] | G.N. Coleman, R. Johnstone, C.P. Yorke, and I.P. Castro. DNS of aircraft wake vortices: the effect of stable stratification on the development of the Crow instability. In V. Armenio, B. Geurts, and J. Frolich, editors, ERCOFTAC Workshop: Direct and Large-Eddy Simulation 7, pages 519-525. Springer, 2008. |
| [5] | R. Corral and J. Jiménez. Fourier/Chebyshev methods for the incompressible Navier-Stokes equations in infinite domains. Journal of Computational Physics, 121:261-270, 1995. |
| [6] | M. Frigo and S.G. Johnson. The design and implementation of FFTW3. Proceedings of the IEEE, 93:216-231, 2005. |
| [7] | J.L. Gustafson. Reevaluating Amdahl's law. Communications of the ACM, 31:532-533, 1988. |
| [8] | N. Hutchins and I. Marusic. Evidence of very long meandering features in the logarithmic region of turbulent boundary layers. Journal of Fluid Mechanics, 579:1-28, 2007. |
| [9] | R. Johnstone, G.N. Coleman, and P.R. Spalart. The resilience of the logarithmic law to pressure gradients: evidence from direct numerical simulation. Journal of Fluid Mechanics, 643:163-175, 2010. |
| [10] | J. Kim, Moin. P., and R. Moser. Turbulence statistics in fully developed channel flow at low reynolds number. Journal of Fluid Mechanics, 177:133-166, 1987. |
| [11] | N. Li. 2DECOMP&FFT - library for 2d pencil decomposition and distributed fast fourier transform. http://www.2decomp.org/. |
| [12] | N. Li and S. Laizet. 2DECOMP&FFT - a highly scalable 2d decomposition library and FFT interface. Cray User Group 2010 conference, Edinburgh, 2010. |
| [13] | L. Snyder. Type architectures, shared memory, and the corollary of modest potential. Annual Review of Computer Science, 1:289-317, 1986. |
| [14] | P.R. Spalart, G.N. Coleman, and R. Johnstone. Retraction: Direct numerical simulation of the Ekman layer: A step in Reynolds number, and cautious support for a log law with a shifted origin [phys fluids 20, 101507, 2008]. Physics of Fluids, 21:109901, 2009. |
| [15] | A. Wray. Vectoral: only by trying new things will we ever get computer languages right. http://merrimac.stanford.edu/brook/vectoral.pdf. |