OpenFOAMスレッド並列化に対する挑戦

TL;DR

• pimpleFoamのDIC-PCGソルバーを、OpenMPを用いたスレッド並列化しました
• DIC-PCGは、スレッド並列化にとって最も難しい解法の1つです
• オープンCAE学会V&V小委員会のチャネル流ベンチマークをIntel KNL上で計測しました
• スレッド並列化することで、非並列版の13.5倍になりました
• 通常のMPI並列版に比べると半分の性能ですが、最も難しい問題においてこれだけの改善であれば、他の問題では同等か優位になるだろうと期待できます
• 引き続き、他のソルバーやベンチマークや環境で研究を進めて、その改善や効果測定を進めていきます

はじめに

動機

現代の機械製品や社会基盤等の設計開発にあたって、計算機援用工学(Computer Aided Engineering)CAEは欠かせないツールです。
そしてその中で、OpenFOAMは、世界で最も使われているCAEツールの1つです。
ですので、特に、OpenFOAMの計算時間を高速化短縮することは、そのような設計開発に必要な解析時間の削減および試行回数の向上のために重要な課題となっています。

一方で、現代の計算機において、OpenFOAMはあまり良い性能を出せていません。
例えば以下の図に示す今野(2017)の結果では、ノード数が増えるにつれて処理速度が低下してしまっています。
この例以外にも、多くの研究で同様の傾向が報告されています。

これは、今のOpenFOAMが現代の（そして将来的に得られるであろう）計算機の性質にあまり合っていないことを意味します。
現代はポストムーアの問題に直面しているとも言われています。
本当にポストムーアの時代が訪れるかはともかくとしても、いずれにせよ、そのハードウェア的な性能がずっと飛躍的に向上されつづける、つまり「ソフトウェアが何もしなくても解析時間は短縮されていく」という状況には期待できません。
よって、現代および未来の計算機でOpenFOAMを効率化・高速化するには、OpenFOAMのソフトウェア的な性質や構造を再検討し改善する必要があります。

背景技術

OpenFOAMの根本的な改善のため、私達は並列手法に着目しました。
今のOpenFOAMは並列化には、MPIのみを用いるフラットMPI並列を採用しています。
一方、昔の計算機と今の計算機の違いを次の表に示します。
この表から、真ん中の列に示すようなOpenFOAMが最初に開発された古い時代に期待されていた計算機（アーキテクチャ）はもはや成立しておらず、今後もそのような前提を期待したフラットMPIを使い続けることは得策ではない事がわかります。

計算機の性質の今昔比較

	古代の計算機	現代の計算機
性質
1ノードあたりのCPUコア数	1コア/少ない	とても多い
ノード数	あまり多くない	とても多い
CPUの演算速度（ネットワーク内通信速度に比べて）	十分に遅い	速い
MPIの文脈では
MPIプロセス数	とても少ない	とても多い
MPIプロセスの管理コスト（演算コストに比べて）	十分に軽い	重い
MPIのプロセス間通信コスト（演算コストに比べて）	十分に軽い	重い

そのような問題を避けるため、OpenFOAMをMPIだけではなく、OpenMPを用いたスレッド並列化もするのが、私達の課題です。
MPIとOpenMPの並列化手法の違いを次の表に示します。
このように、ノード内の並列であれば、スレッドを用いた共有メモリ型並列を使うことで性能向上が期待できます。

並列化手法の比較

	今のOpenFOAM	私達の方法
利用するフレームワーク	MPI	OpenMP
実装方法	プロセス	スレッド（軽量プロセス）
主なデータ通信方法	ソケット	共有メモリ
利用できる主な環境	ノード内も外も可	同じノード内のコア間
コストの面では
管理コスト	重い	軽い
通信コスト	重い	軽い

ということで、私達の目標は

• OpenFOAMに、OpenMPを用いた並列化も実装し、
• ノード内はスレッド並列、ノード間はプロセス並列となるハイブリッド並列を実現させ、
• その改善効果を評価し、
• その結果や実装を多くの人に共有する

ことになります。

具体的な手法

具体的には、これまでで、私達は以下のようなことに取り組みました。

• チャネル流れのベンチマーク（channelReTau110問題）を解くため、
• 疑似SIMPLE法（pimpleFoamソルバー）の
• 対角不完全コレスキー分解前処理付共役勾配法（DIC-PCG線形ソルバー）を用いて、
• メニーコアCPU (Intel Xeon Phi Knights Landing)上で
• 1ノードでOpenMPを使って並列化・高速化

OpenFOAMで解ける問題は非常に多様なのですが、今回は、オープンCAE学会のV&V小委員会が提供しているもので、先述の今野(2017)が使用しているものと同じものをまずは対象としました。
pimpleFoamは、このベンチマークで使われている流体ソルバーです。
今回はpimpleFoam内の処理の内、特に線形方程式の求解処理に注力しました。これは、多くの（陰解法型）流体計算において線形方程式の求解に最も時間がかかっていることが経験的に知られているからです（後述の通りpimpleFoamでの実測結果もそのような結果となっています）。
線形ソルバーの選択にも余地があるのですが、「まずは、スレッド並列化の効果が最悪でもどれぐらい得られそうか」を知るために、スレッド並列にとって最も困難な問題の1つであるDIC-PCGにしました。
計算機については、先述の通りフラットMPIではコア数の多さも問題になるため、コア数が多いKNLを用いました（今現時点においてKNLは少し時代遅れ感がありますが、将来的に他のCPUも、KNLと同様な方向にコア数が増えていくことが予想されるため、その状態を予見してあえて使います）。
今回の実装はハイブリッド並列になるので、もちろん複数ノード実行もできるのですが、問題を単純化するためにノード数は1つで実行・評価します。

この状態で期待されるのは、フラットMPI版に比べて同等か最低でも少し遅いぐらいです。そうすると、複数ノードにした時に、フラットMPIより速くなることが期待できます。

以下では、実際にDIC-PCGをどのようにスレッド並列化したかについて示します。

疎行列ベクトル積

通常の行列とベクトルの積を思い出してもらえば分かる通り、この演算は基本的に行ごとに独立のはずです。
しかしながら、行列が疎（多くの要素の値が0）である「疎行列」においては、行列のデータ量・演算量を減らすために工夫された格納形式を持っています。
pimpleFoamの使うDIC-PCGでは、lduMatrixという形式を使っています。

この形式はpimpleFoamに限らず多くの流体ソルバーでも使われているもので、例えば
\[\begin{pmatrix}
1 & & 5 & 6 \\
& 2 & & \\
8 & & 3 & 7 \\
9 & & 10 & 4 \\
\end{pmatrix} \]
という行列があった時、lduMatrixは次の表のようなフィールドを使って表現されています。

lduMatrixの例

フィールド名	値	意味
diag	1,2,3,4	対角要素の値
upper	5,6,7	上三角行列の値
lower	8,9,10	下三角行列の値
upperAdder	2,3,3	上三角行列の列番号（下三角行列の行番号）
lowerAdder	0,0,2	下三角行列の列番号（上三角行列の行番号）

このように、lduMatrixでは

• 下三角行列L（列優先）
• 対行列D
• 上三角行列U（行優先）

という3つに分けて値を保持しており、この形式は一般的にCOO形式と呼ばれるものです。
そして、この形式において、下三角行列Lが列優先であることが、疎行列ベクトル積の並列化を妨げています。

for (label face=0; face<nFaces; face++)
{
	ApsiPtr[uPtr[face]] += lowerPtr[face]*psiPtr[lPtr[face]];
	ApsiPtr[lPtr[face]] += upperPtr[face]*psiPtr[uPtr[face]];
}

lduMatrixのAmulにおける変数の具体例（先述の行列を例に）

face	0	1	2
uPtr	2	3	3
lPtr	0	0	2

上のコードと表は、実際の疎行列ベクトル積のコードと、各反復変数faceの時のuPtr(upperAdder)とlPtr(lowerAdder)の具体例です。
コードを見ると分かる通り、ApsiPtrの書き込み先はlPtrとuPtrの値によって決まりますが、この表で例えばlPtrがface=0とface=1で同じになっているように、別の反復変数faceであっても、書き込み先が同じになることがあります。

このような反復を単純に#pragma omp parallel forと入れてOpenMPでスレッド並列化してしまうと、別々のスレッドが同時に同じ変数に書き込むことになり、データ競合というバグになります。

その問題を解決するため、疎行列の格納形式をCSR形式に変更しました。
CSR形式は、疎行列の値を、行方向に「値」「列番号」と「その行の先頭位置」で格納するもので、先の行列の例では以下のようになります。

CSR行列の例

フィールド名	値	意味
data	1,5,6, 2, 8,3,7, 9,10,4	非ゼロ要素の値
column	0,2,3, 1, 0,2,3, 0,2,3	各要素の列番号
offset	0, 3, 4, 7, 10	各行の先頭要素の番号

CSR形式を使った疎行列ベクトル積は以下のようになります。このように、各行で書き込み先が独立になるため（反復変数iに対してy[i]のみ）、データ競合は発生せず並列処理が可能になります。

for (label i = 0; i < n; i++)
{
	double y_i = 0.0;
	for (label index = offset[i]; index < offset[i + 1]; index++) {
		y_i += data[index] * x[column[index]];
	} 
	y[i] = y_i;
}

DIC前処理における前進後退代入

DIC前処理は、名前の通り、不完全コレスキー分解、もっと言えばLU分解の一種であり、

1. 最初に元の行列AをA=LUという下三角と上三角行列の積となるようなLとUに分解し
2. Lを使って前進代入
3. Uを使って後退代入

をする処理です。この前進・後退代入が並列化の時に問題になります。

for (label face=0; face<nFaces; face++)
{
	wAPtr[uPtr[face]] -= rDPtr[uPtr[face]]*upperPtr[face]*wAPtr[lPtr[face]];
}

上のコードは、実際の前進代入の箇所のコードですが、ご覧の通り、左辺と右辺の両方にwAPtrが出現しています。
ですので、例えば例題の行列だと、face=2の時にlPtr=2なのですが、wAPtr[2]が更新されるのは、face=0の時なのでface=2を計算する時にface=0を先に計算していなければいけません。
もしこれをそのまま並列化してしまうと、スレッドの実行順序によってfaceのどちらが実行されるか分からないため、実行毎に結果が変わってしまう競合状態が発生します（厳密には同時に読み書きが発生するためデータ競合にもなりうるのですが、それはwAを適切にアトミック変数にすることで回避が可能です。しかしアトミック変数にしても競合状態は回避できません）。

これを回避するには、計算順序を変える順序並び替えと呼ばれる特別な技法が必要です。
今回、私達はその並び替えにCuthill-McKee法を用いました。
Cuthill-McKee法による並び替えでは、下図のように依存関係を保ったまま並列化できるため、他の順序並び替え手法（例えば多色順序付けなど）に比べて収束性の劣化が抑えられることが知られています。

一方で、この依存関係の抽出処理に時間がかかる上、並列数は独立した行（セル）の数に上限があるため並列化効率はあまり上がりません。
これが、DIC前処理がスレッド並列化にとっても最も難しい問題の1つとなっている理由です。

Cuthill-McKee法を用いた詳細については、中島研吾(2013)などが詳しいのでそちらをご参照ください。

処理時間の計測結果

以上のようなOpenMPを用いたスレッド並列化版の改善を評価するため、以下の条件で処理時間を計測しました。

ベンチマーク問題

オープンCAE学会V&V小委員会の”channelReTau110″

ソルバー

DIC-PCG線形ソルバーを選択したpimpleFoam

元となったOpenFOAMバージョン

Foundation版（コミットハッシュ16b559c1）

計算機

Intel Ninja Developer Platform（Intel Xeon Phi 7210; 256論理コア、DDR4）

計測方法

最初の5ステップの平均値。後述の詳細プロファイルは各関数の先頭と最後に手動でタイマーを仕込んで算出した

非並列版での計測

まず、並列化されていない版での状態を知るために、何も変更していない版で実行時間の内訳を計測しました。結果は上図のようになり、

• 先述の通り、pimpleFoam全体では線形方程式の求解処理がほとんどを占めていること
• PCGの中でも、AmulとDICが最も時間がかかっていること

を確認しました

スレッド並列化の効果計測

上図は、先述のスレッド並列化技法を順番に適用していった結果です。

1. 一番左の結果は元となる非並列版の結果で、先の図の一番右と同じです。
2. 疎行列格納形式をCSRに変えた結果です。結果、DIC前処理が遅くなってしまいました。これは、キャッシュヒット率の低下によるものと考えられます。
3. そこで、疎行列格納形式を少し変えました。その方式は、CSRの形を保持したまま、値をldu形式のように上三角・下三角・対角の3つに分けて保持するもので、ldu-CSR形式と呼ぶことにします。これにより性能劣化は解消されもとに戻りました。
4. 準備は整ったのでここから並列化に入ります。まずは先述の通り最も時間がかかっている行列演算を並列化しました。高速化率はAmulがx12.8、DICがx3.5で、全体として約2倍になりました。ここでCuthill-McKeeの計算時間が追加され、全体の2割を占めることになりました。
5. 行列演算を高速化すると、次はベクトル演算の時間が目立つようになるので、こちらも高速化します。WAXPBYは50.3倍、sumMagは13倍、sumProgは13.8倍で、全体で元のバージョンから3.4倍になりました。
6. 以上の結果はOpenMPの既定の設定のままであり、スレッド数が256、つまり1物理コアに対して4スレッドのハイパースレッディングになっていました。計測した結果、1コア1スレッド、つまりハイパースレッディングせず64スレッドが最も高速で、元のバージョンに比べて4.8倍になりました。

ということで、最終的に全部合わせると4.8倍なのですが、ここで、本来はCuthill-McKee法の時間については無視することができます。
現在の実装では、毎回の時間反復のたびにCuthill-McKeeをしていますが、これは無駄で、メッシュが再構築された時（つまり今回のベンチマークであれば初回）に1回実行すれば済みます。
かつ、実際の計算では5ステップではなくもっと長時間の計算を実行しますから、相対的に少ない回数しか実行しないCuthill-McKee法は無視してよいはず、ということです。

ということで、Cuthill-McKee法が無視できたと仮定して、最終的には13.5倍の加速率を得ることが出来ました。

各関数の処理時間と加速率

関数名	非並列版[s]	OpenMP版[s]	加速率[s]
Amul	60.0	3.0	x19.8
DIC::precondition	80.1	7.9	x10.1
WAXPBY	21.1	0.4	x53.2
sumMag	6.1	0.1	x44.6
sumProd	12.3	0.3	x43.6
Cuthill-McKee	0.0	24.0	—
(other)	1.1	1.6	x0.7
total	180.8	37.4	x4.8
total excluding Cuthill-McKee	180.8	13.4	x13.5

フラットMPIとの比較

スレッド並列化によって並列化する前に比べて13.5倍を得ましたが、では元々のフラットMPIではどうだったのでしょうか。

フラットMPIとの比較結果を上図に載せます。
ご覧の通り、今回実装したOpenMPによるスレッド並列化版は、Cuthill-McKee法を抜いても、フラットMPIに比べて半分の処理性能（倍の処理時間）となってしまいました。

この原因については2つ挙げることができます。1つ目の理由は、下の表のように、DIC前処理だけが、他に比べてあまり高速化されていなかったからと考えられます。

DICの加速率があまり良くない原因は2つあります。1つは、先述の通り、理論的に、どうやっても並列数が頭打ちになってしまうという元から分かっていた問題。
もう1つは実装上の問題で、現在、入出力ベクトルの並び替え処理を毎回やっているのですが、これは原理的には最初に1回だけやればいいので、そのように変更すると改善されると考えられます。

また、OpenMPがあまり速くないもう1つの理由は、今回用いたベンチマークの性質にあります。
本来であれば、MPIプロセス数を増やすと、領域分割数が増えるため、収束性が悪化し反復回数が増えるはずです。
一方、OpenMPによるスレッド並列化であれば、領域分割法を使わないため、反復回数は元のままです。
これによりOpenMPの方が速くなることが期待されたのですが、今回用いたベンチマーク（チャネル流）は、規則正しい正規格子であり、流れ場も一様という、非常に行儀の良い問題であったため、領域分割による性能悪化が現れなかったと考えられます。

このような理由により、処理速度は残念ながら半分程度になってしまいましたが、逆に言えば、これだけ条件の悪い条件で「半分」で済んでいるという結果とも言えます。
ですので、今回の結果が示すのは「OpenMPによるOpenFOAMのスレッド並列化は、十分に改善の可能性が期待できる」ということだと私達は結論づけました。

まとめと今後

ということで、ここまでの成果をまとめると、OpenFOAMをスレッド並列化し

• 非並列版に比べて、13.5倍の高速化
• フラットMPI版に比べても、半分程度に迫る高速化

という成果を、「正規格子で一様場という単純な問題」で「DICという最も難しい手法の1つ」を使って得ることが出来ました。

冒頭で述べた通り、今やっている取り組みはまだまだ途中であり、今回の内容は1つの得られた結果報告に過ぎません。今後は

• Cuthill-McKee法自体を並列化し、再メッシュ切りされた時にも高速化の恩恵が受けられるようにしたい
• DIC前処理はもう少し高速化したい（例えば、入出力ベクトルの並び替えの回数削減）
• 他の前処理やソルバーの並列化と高速化性能測定も実施し、OpenMP並列のその他の問題点を見つけたい
• もっと複雑なメッシュ条件を用いたベンチマーク問題を解いて、領域分割がないことによる効果によって、OpenMP版の方が速くなるはずという期待が正しいか調査したい
• 複数ノードの計算機、特にスーパーコンピューターにおける性能調査で、MPIプロセスが減ることによるコスト削減の効果も調べたい

と、まだまだやり残したことがたくさんあり、引き続き、これらに取り組んでいきたいと思っています。また結果を得られ次第、本ブログ等でご報告いたしますので、お楽しみに！

また、今回実装したコードについては、もう間もなく公開する予定です（現在コード等の整理中）。コード公開できましたら、そちらもお知らせします！

最後に

今回の成果は私達のお客様から具体的な発注があって納品した成果ではなく、社内有志が進めていたプロジェクトの成果になります。
となると、「なんでフィックスターズがOpenFOAMをやってるの？」と疑問に思われるかもしれません。
なぜなら、フィックスターズは流体解析を仕事にしているOpenFOAMのユーザーではなく、OpenFOAMやその他の流体ソルバー・ツールを売っている会社でもないからです。

端的に言えば、これは「ソフトウェアの高速化」が本業であるフィックスターズにとっての、高速化実例のひとつとして取り上げたに過ぎません。
そのような実例は、1つは社内のエンジニアに対して技術研鑽の課題として、1つは（潜在的な）お客様に対しての具体的なサンプル提供となります。
そういった（個人的な）動機もあり、今回はあえて難しい課題から挑戦してみた結果、となります。

ということで、フィックスターズは、こんな難しい課題を解いて楽しむ文化のある会社です。
そんな会社で働くことに興味がちょっとでも沸いたなら、会社説明会などにぜひお越しください。

また、OpenFOAMや流体計算はもちろん、それ以外にも難しいものから簡単なものまで「ソフトウェアの処理時間を短くしたい」というご要望があれば、ぜひご相談ください。
もちろん、OpenFOAMの高速化そのものに興味がある方からのお問い合わせも大歓迎です！

ライセンス

本記事は、発表資料と同様に、OpenFOAMのコードなどを含んでいるため、OpenFOAMと同様にGPL v3で利用許諾されています。
利用に際しては、以下の表記を用いて適切にご利用ください。

GNU GENERAL PUBLIC LICENSE

This article contains OpenFOAM source code.

OpenFOAM is Copyright (C) 2011-2017 OpenFOAM Foundation
Contact: OpenFOAM Foundation, http://openfoam.org/contact

Otherwise noted, other texts and images in this article are Copyright (C) 2019 Fixstars Corporation.

You may use, distribute and copy this article under the terms of GNU General Public License version 3, which is displayed at https://www.gnu.org/licenses/gpl-3.0.html, or (at your option) any later version.