RIAM-COMPACT HPC 版の開発
九州大学 応用力学研究所

内田 孝紀

九州大学 情報基盤研究開発センター 小野 謙二
1. 開発方針
風車の設置場所を選定するための風況予測シミュレーションは，事業の成功を左右する判断材料となるため，主に下記の観点で高精度な予測結果が要求される．

1．複数の風向での評価（16−32 方位）2．10 分間の実時間平均値
3．様々な風況での評価（乱流強度や変動風など）

次世代の風況予測では，従来のシミュレーションのように数点の評価ケースから判断を下すのではなく，これらの点を考慮したシミュレーションを行い，多数の結果から総合的な評価を行う必要がある．そのためには，利用する計算機の性能を引き出せるシミュレータを構築し，多数の計算を短時間で実施することが求められる． 利用する計算機は，想定するユーザ層を考慮して，富岳を始めとして Intel のクラスタ，ベクトル機など多様な計算機を検討しなければならない．そのため，異なるアーキテクチャでも高速な計算が可能な計算方法が望ましい．
基礎方程式として，非圧縮の Navier-Stokes 方程式を用いるため，圧力の Poisson 方程式の求解が最も時間のかかる部分となる．したがって，大規模疎行列を係数行列にもつ反復解法の高速化が最重要課題の一つとなる．
2. 多風向同時計算の検討
プロジェクト開始時には，富岳のアーキテクチャの方向として，低Ｂ/F のメニーコアであることが想定されたため，低Ｂ/F アーキテクチャで高速な実行が可能な反復法の開発が主な目的の一つであった．ループ内の演算密度を増加するアイデアとして，疎行列ベクトル積（Ax=b）を行列行列積の形式(AX=B)に変換することを検討した．つまり，同一の係数行列に対して異なる右辺ベクトルを用意し，対応する複数の解ベクトル X を同時に求める方法である．異なる右辺ベクトルの具体的な例は，異なる風向の計算が考えられる．この方法は，利用する係数行列をメモリからキャッシュに一度ロードすれば，その値を複数の解ベクトルの求解に流用できるため，要求 B/F を減少させることができ，低 B/F アーキテクチャでの高速化が期待できる．
まず，異なる風向の計算を同時に実行する場合，流れの状況によって計算時間が異なる可能性がある．つまり，流れの変動の度合いにより圧力場を求める場合の反復回数に違いが生じる．速度場の影響は連立一次方程式の右辺項に含まれるが，速度場に変化があると右辺項の値が変化するため，それに応じて収束性が変化することがわかっている．風車の場合は，風向や配置により風車間の相1互干渉の程度が変わるため，前述のような状況を引き起こす可能性がある．これを確認するために，
簡単な実験を行った．
実験は，二次元の 2 角柱を流入風向に対して直列に配置し，角柱間の距離を変化させたシミュレーションを行った．シミュレータには，非圧縮性流体シミュレータである FFV-C を利用した．FFV-Cは直交系格子を用いた LES シミュレータであり，開発する RIAM-COMAPCT と類似の特性を示す．図 1に各中間の中心距離ｘを x=1.2〜4.5 まで変化させた場合の流れの様子を計算し，各ケースの圧力反復回数の履歴を示す．横軸はステップ数であり，縦軸は反復回数を示す．配置により，流れの状況が変わるため，圧力の反復回数の履歴が異なることが確認できる．

図 1 異なる配置に対する流体シミュレーションの圧力反復回数の履歴
図 2 には，最も反復回数の少ない x=1.5 のケースを基準に 30000 ステップまでの反復回数の総数を比較した．これより，今回の単純なケースでも，35％程度の反復回数の差が生じることがわかる．

2

図2

圧力反復回数の比較

したがって，行列積行列積の形式で要求 B/F を低減し，高速化を図ることは可能であるが，この方針を複数の風向として実装する方法は全体として効率化することが保証できないという結果となった．そこで，別の方法として，オーソドックスに 1 つケースの計算時間を短縮することを考える．1
ケースの計算を短時間化できれば，多風向の計算を同時に解かずとも同等の効果を得ることができる．
3. 低 B/F 疎行列反復解法
反復解法部分の要求 B/F を低減するために，高 B/F となる疎行列反復法に部分的に直接法を組み合わせた方法を提案した．提案方法は，L1 キャッシュの利用効率が非常に高く，メニーコアアーキテクチャにおけるスレッドスケーラビリティが非常に高い特性を持つ．また，提案手法は，メニー
コア，マルチコアアーキテクチャだけでなく，ベクトル計算機においても非常に高い性能を達成できることが確認できた．
圧力の Poisson 方程式を離散化すると，疎行列を係数行列にもつ連立一次方程式が得られる．式（1）では物理量φ，格子幅 h，ψを Poisson 方程式のソース項としている．

(1)
式(1)に対して，古典的な SOR 反復法を適用すると，

3

(2)ここでハット記号(^)は k 方向に関する三重対角行列を直接法で解いた後のテンポラリの解ベクトルである．その後，式(3)の SOR の反復過程を適用する．

(3)
式(2)の三重対角行列は，最も効率的な LU 分解（トーマスアルゴリズム）を用いて解くことが多い．
ただし，LU 分解は逐次アルゴリズムであるため，並列化も SIMD 化もできない．また，効率的な省演算解法で高い B/F アルゴリズムであるため，高い演算性能が得られない．そこで，LU 分解を並列化するため，Parallel Cyclic Reduction (PCR)を採用する．
3.1 PCR の概要LU 分解は前進消去と後退代入の 2 パスのアルゴリズムである．各パスで，依存関係があるために並列化ができない．PCR 法はその依存関係がなくなるように，式変形を行う．1 回の式変形操作を縮約という．1 回縮約すると，最初の N 元の連立一次方程式は，2 つの N/2 元の連立一次方程式となる．
更に再度縮約操作を行うと，4 つの N/4 元連立一次方程式が得られる．これを繰り返すと，最終的にN つの単一方程式が得られ，直ちに解が得られる．
連立一次方程式 Au=d を考える．

(4)
この行列の i 行を考えると，4(5)式(5)は点 i について，i±1 の点が依存関係をもつ．i±1 行の式に行列の基本操作を施し，i±1 点との依存関係を断ち切るようにすると，(6)が得られる．式(6)は両隣の点と直接関係がないことがわかる．この処理を縮約処理とよぶ．なお，この縮約操作により，係数の値が変化していくことに注意する．再度，縮約処理を施すと次式を得る．

(7)
縮約処理は全ての点との依存関係がなくなるまで繰り返す．処理の流れを図 3 に示す．

図3

PCR アルゴリズムの縮約処理．色は異なる方程式の組を示す．

解法の詳細については，[1] を参考にされたい．
3.2 評価機
SLOR-PCR 法を評価する計算機として，マルチコア Intel Skylake-SP, メニーコア SGI UV300, ベ
クトル NEC SX-Aurora の 3 機種を選定した．主要諸元を表 1 に示す．

5

表 1 評価機の諸元

特に，UV300 は富岳のアーキテクチャと類似した構成となっている．つまり，メニーコア，cc-NUMAアーキテクチャ，リングバス構造などである．UV300 での特性を考察することにより，富岳での動作の傾向を推測することができるだろう．

富岳

UV300

http://pr.fujitsu.com/jp/news/2018/08/22-1.html

図4

富岳と UV300 のアーキテクチャの比較

3.3 評価結果
まず，図 5 に Skylake-SP における SLOR-PCR 法の演算性能を Jacobi 反復法，Red-Black SOR 法と比較して示す．SLOR-PCR 法は３６コアまでスケーラブルなスレッド性能を示している．一方で，古典的な反復解法では演算性能は低い．両者の差は約６倍程度あることがわかる．
6

図5

Skylake-SP における反復性能の比較

次に，演算性能の違いの原因を探るため，キャッシュの利用率を調べた．図 6 に L1 キャッシュとL2 キャッシュのヒット率を示す．SLOR-PCR の L1 ヒット率は 98%知恵度であるのに対し，Red-BlackSOR 法は 76%程度である．一方，L2 ヒット率については SLOR-PCR 法は 7%程度でほぼ一定であり，L1キャッシュが効率的に効いていることがわかる．Red-Black SOR 法は L2 のヒット率は 25~60%とばらつきがあり，L1 のヒット率が低いため，L2 まで見ていることがわかる．このキャッシュの挙動により，SLOR-PCR 法は優れた L1 上データの再利用性があることがわかる．

図 6 Skylake-SP におけるキャッシュヒット率7図 7 に示すように，アーキテクチャの異なる SX-Aurora についても，同様の挙動を示すことがわかった．

(a) 演算性能(b) キャッシュヒット率

図 7 SX-Aurora における演算性能とキャッシュ挙動
図 8 には UV300 の演算性能を示す．UV300 は 4 つのユニットが NumaLink7 で接続されており，各ユニットは 88 コアをもつ．全 352 コアで SLOR-PCR 法は他の 2 つの反復解法に比べて高い演算性能を達成している．特に，1 ユニット内では 88 コア時に 80 倍以上のスピードアップを達成しており，メニーコアの性能を十分に引き出せていることが確認できる．

図8UV300 の演算性能の比較

83.4 計算時間の比較
これまでは，計算速度に対する評価を行ってきたが，図 9 には実際の計算時間を評価する．図 9は右側が解きにくい問題であり，左側が解きやすい問題を示す．Jacobi 法やRed-Black SOR 法は問題が解きにくくなると急速に反復回数が増加し計算時間が増えるが，SLOR-PCR 法の反復回数は増加しにくい．これは，直接法における求解の効果であると考えられる．

図9

異なる右辺項に対する収束特性の時間比較

4. RIAM-COMPACT HPC 版による Alpha Ventus 洋上ウィンドファーム解析
次に，今回開発した数値風況予測シミュレータ RIAM-COMPACT HPC 版を用いたベンチマークテストに関する研究成果を報告する。
本研究にてターゲットにしたサイトは図 10 に示す Alpha Ventus 洋上ウィンドファームである。ここは，ドイツ初の本格的なオフショア・ウインドファームであり，5MW クラスの風車が合計 12 基設置されている。

9

図 10 Alpha Ventus 洋上ウィンドファーム
図 11 には，本研究で対象とした風車の諸元を示している。これらの数値は空力弾性シミュレーション
コード FAST にて得られた結果である。図 12 には，本プログラムの検証に用いた工学モデル(石原モデル)の結果を示す。

図 11 NREL 5MW wind turbine の仕様

10図 12 工学モデル(石原モデル)の概要
最初に，風車単機の計算を行い，予測精度の検証を実施した。その結果を図 13 に示す。Near wake 領域および Far wake 領域ともに概ね良好な一致を得ることに成功した。

図 13 ベンチマーク計算の結果，工学モデル(石原モデル)との比較

11

図 14 ベンチマーク計算の結果(風車 3 基の場合)，工学モデル(石原モデル)との比較
次に，風車を 3 基設定した場合の検討を行った。得られた結果を図 14 に示す。特に，下流 3 台目の風車の下流域に注目する。風車の相互干渉の影響を強く受け，工学モデルと有意な差異が生じていることが明確に観察される。この結果は，既存の工学モデルは複数台の風車ウエイクの相互干渉の再現には問題あることが分かる。
実スケールの Alpha Ventus 洋上ウィンドファームを対象した計算を実施する前に，最終的な計算コードの確認を，以下に示すミニ Alpha Ventus 洋上ウィンドファームを作成して行った(図 15 を参照)。今回のミニ Alpha Ventus 洋上ウィンドファームでは，風車のスケールおよび台数は，実際の Alpha Ventus
洋上ウィンドファームと同じである。異なるのは風車の離隔距離である。ミニ Alpha Ventus 洋上ウィンドファームでは風車の離隔距離を，風車直径 D の 3 倍と設定した。計算格子は主流方向，スパン方向，鉛直方向にすべて同一とし，2m に設定した。全格子点数は 8,800 万点である。風車ブレードの回転にはアクチュエータラインモデルを採用している。風車の周速比は 7.0 である。流入風の傾きは N=10 とした。
図 16 および図 17 には，ミニ Alpha Ventus 洋上ウィンドファームを対象とした計算結果例を示す。特に図 17 では，計算プログラムの有効性および妥当性を検証するために，任意の風車配置によるテストを実施した。両者のケースともに計算は安定に実行された。流れ場の可視化結果を吟味すると，複数の風車下流に形成されている風車ウエイクが相互に干渉している様子が確認される。形成された流れ場は，強非線形な流動現象であり，我々が採用した乱流モデル LES によるアプローチの有効性が示された。

12

図 15 ミニ Alpha Ventus 洋上ウィンドファーム

図 16 ミニ Alpha Ventus 洋上ウィンドファームを対象とした計算結果例①13

図 17 ミニ Alpha Ventus 洋上ウィンドファームを対象とした計算結果例②

最後に，実スケールの Alpha Ventus 洋上ウィンドファーム(風向は西風)を対象にした計算結果を示す。
空間解像度は各方向ともに 5m である。総格子点数は 821(x)*721(y)*231(z)である。約一億四千万点である。図 18 には，風車ハブ高さ 92m における風速分布(瞬間場)を示す。各風車から風車ウエイクが形成され，下流側ではそれらが複雑に相互干渉している様子が明確に観察される。
特に下流 3 列目の風車群(No.3,6, 9, 12)は上流側に位置する 2 台の風車ウエイクの影響を受けている。図 19 には，風車ハブ高さにおける風速情報とそのグラフを示す。各風車位置にて特徴的な傾向が再現されている。

14

図 18 Alpha Ventus 洋上ウィンドファームを対象とした計算結果例15

図 19 Alpha Ventus 洋上ウィンドファームを対象とした計算結果、風車ハブ高さの風速情報

5. まとめ以上に述べたように，直接法と反復解法のハイブリッドである SLOR-PCR 法は L1 データの再利用性が高く，メニーコアやベクトルアーキテクチャに対して高性能を達成できる反復解法であり，その性能は古典的な反復法に比べて 6 倍程度であることがわかった．
さらに，直交不等間隔格子を採用した場合に生じるメトリクス計算部分を考慮すると，さらにループ内の演算密度を高めることができ，コードの要求 B/F を下げることができる．

16

今後，RIAM-COMPACT HPC 版は，富岳のみならず，多様な計算機での利用を念頭に整備開発を継続していく．
また同時に，実スケール Alpha Ventus 洋上ウィンドファームの解析については，他の風向などの計算を行い，総合的に RIAM-COMPACT HPC 版の性能を検証していく予定である。

参考文献[1] Ono, K., Kato, T., Ohshima, S., and Nanri, T. : Scalable Direct-Iterative Hybrid Solver forSparse Matrices on Multi-Core and Vector Architectures, Proceedings of the International
Conference on High Performance Computing in Asia-Pacific Region, pp. 11–21,
ACM HPCAsia2020, Fukuoka, Japan, 2020.

17