J. Plasma Fusion Res. Vol.90, No.12 (2014)７
５５‐７６３

解説

GPU コンピューティングによる大規模シミュレーション
Large-Scale Numerical Simulations Based on GPU Computing
青木尊之
AOKI Takayuki
東京工業大学 学術国際情報センター（原稿受付：２
０
１
４年７月１
６日）

GPU は高い演算性能とメモリアクセス性能をもつアクセラレータとして多くのスパコンに使われている．さらに低消費電力であり利点が多いが，既存のプログラムを書き換える必要がある．内部に２０００個を超える演算ユ
ニットをもつ GPU のアーキテクチャの説明，並列計算を前提とする CUDA プログラミングの概略を述べる．また，通常のプログラムに指示行を追加するだけで GPU を利用できるようにする OpenACC についても簡単に触れる．著者が GPU スパコン TSUBAME を利用して行った大規模シミュレーションの例として，水平解像度 500 mの気象計算，1 m 格子を用いた 10 km 四方の東京都心部の超高解像度気流計算，半陰解法による気液二相流計
算，３３００億格子を用いたフェーズフィールド法による合金の樹枝状凝固計算，離散要素法による粉体シミュレー
ションについて解説する．最後にプラズマ・シミュレーションにGPUコンピューティングを適用する見通しについて示す．
Keywords:
GPU, supercomputer, TSUBAME, CUDA, weather prediction, turbulent LES, gas-liquid two-phase flow,
phase-field model, dendritic solidification, DEM, dynamic load balance

１．GPU について

れ，GPU の性能向上に合わせた頻繁なバージョン・アップ

１．
１ GPU の汎用計算への適用とともに機能もますます向上している．
１．
２ GPU スパコンと Top500 ランキング

パソコンのグラフィクス・ボードや CPU，最近ではスマートフォンやタブレットの中にも画像処理プロセッサ

最新の GPU は単体で 1.7 TFlops（倍精度演算），単精度

GPU（Graphics Processing Unit）が搭載されている．GPU

演算では 5.3 TFlops を超えるようなピーク演算性能を持は元々画像を表示するために開発され，その後，より高速

つ．また，スパコン専用に開発されたプロセッサと違い，に，美しく，高精細に画像を表示するために描画機能が飛GPU はパソコンやモバイル IT の巨大な市場で展開する製

躍的に向上してきた．そして，画像表示だけではなく汎用品であるためにコストが低い．

計算に使う「GPUコンピューティング」（GPGPU（General-

２００８年に東京工業大学・学術国際情報センターが６
８０個

Purpose computing on GPUs）と同義）の取り組みが２０００年の GPU をスパコン TSUBAME 1.2 に導入し，世界から大き頃から始まった．GPU はグラフィクス処理を目的としていな注目を集めた．NVIDIA は２
０１０年に倍精度浮動小数点演るため，汎用計算に対する機能が限定されている反面，同算性能の向上や ECC メモリに対応した Fermi コアの GPU

世代の CPU に対して１０倍以上高い演算処理能力をもつたをリリースし，GPU をスパコン分野で利用できる条件がめ，魅力的な演算デバイスとして次第に認識されるように整った．２
０１０年１１月のスパコン Top500 のランキングでなっていった．２００６年には NVIDIA 社が自社の GPU に対は，１位と３位に中国の GPU スパコンが入り，４位にはし，GPUコンピューティングの統合開発環境としてCUDA

NVIDIA Tesla M2050 を４２２４個搭載した東京工業大学・学（Compute Unified Device Architecture）を リ リ ー ス した術国際情報センターのスパコン TSUBAME 2.0（総合演算［１］．それまでは Cg（C for Graphics）言語や HLSL（High性能 2.4 PFlops）がランクインし，GPU スパコンが上位をLevel Shader Language）により画像処理の機能を汎用計独占した（図１）．

算に置き換える必要があったが，CUDA の登場により通常

TSUBAME 2.0 の 設 置 面 積 は３３５平 方 メ ー ト ル と 少 なの C 言語の拡張としてプログラミングすることができるよ

く，演算ノードは直前まで事務室だった２つの部屋に収うになり，GPU を汎用計算のアクセラレータとして利用すまっている．２０１２年１１月には，Oak Ridge 国立研究所の Op-

る GPU コンピューティングが一気に広まり始めた［２］．まteron CPU ベースのスパコン Jaguar（２００９年１１月〜２０１０年

た，２００９年 に は CUDA の FORTRAN 版 も リ リ ー ス さ

６月 ま で 世 界１位）に１８６８８個 の Kepler コ ア の GPU

Global Scientific Information and Computing Center, Tokyo Institute of Technology, TOKYO 152-8550, Japane-mail: taoki@gsic.titech.ac.jp

７５５!2014 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.90, No.12 December 2014

試 作 機 TSUBAME-KFC が 4.4 GFlops/W で１位 と な り，
TSUBAME2.5が2.9 GFlops/Wで８位，「京コンピュータ」
は 0.8 GFlops/で１
２５位となっている．Green500 のトップ１
０のスパコンは全て NVIDIA Tesla K20X を導入した GPU
スパコンであり，スパコンの省電力化にGPUが大きく貢献
していることが分かる．
１．
４ GPU のアーキテクチャ
GPU のアーキテクチャは汎用 CPU と大きく異な る．
図１

NVIDIA の Tesla K20X を例 に 説 明 す る（図２）．Kepler

東京工業大学・学術国際情報センターの TSUBAME.

コアの GPU では，１９２個の演算ユニット（CUDA コア）が（NVIDIA Tesla K20X）を導入した Titan が世界１位となっ１つのストリーミング・マルチプロセッサ（SMX と呼ばた．東京工業大学・学術国際情報センターでは，２０
１３年９れている）を構成している．Tesla K20X には１４個の SMX

月に４２２４個の GPUを全て Tesla K20Xに交換し，TSUBAMEがあり，合計で２６８８個の CUDA コアがある．１つの SMX

2.5 として総合演算性能を 5.7 PFlops に向上させている．

の中には 64 k 個の 32 bit レジスタファイルがあり，非常にスパコンに導入される GPU は通常の CPU のマザーボー

高速な L1 キャッシュと共有メモリ（合計で 64 kB），sin

ド上の PCI Express バスに装着されるビデオカード（グラや cos などの数学関数を高速に計算する SFU（Super Func-

フィクス・カード）上に搭載され，アクセラレータと呼ば

tion Unit）などがある．１つの SMX 内の共有メモリは１９２れることも多い．Intel は２０１２年に MIC アーキテクチャの個の CUDA コアからアクセスできるが，他の SMX 内の

Xeon Phi（コードネーム

Knight Corner）を市場にリリー

CUDA コアからはアクセスできない．K２０X を搭載したグスしている．Xeon Phi 5120P はピーク浮動小数点演算性能

ラフィクス・ボード上には 6 GB のビデオメモリ（GDDR5）

1.0 TFlops（倍精度演 算）と メ モ リ バ ン ド 幅 320 GB/sec

があり，GPU からビデオメモリへは 250 GB/sec というメのカタログ性能をもち，GPU とほぼ同等の性能を持つ．描

インメモリに比べて数倍以上高速なメモリアクセスが可能

画機能はもたないが，Xeon Phi もマザーボードの PCI Ex-

である．ビデオメモリには２６８８個の CUDA コアから共有でpress バスに装着して計算の高速化を目的とするアクセラ

メモリアクセスが可能である．通常はCPUからビデオメモ

レータであり，機能を制限した通常のCPUコアを数１
０個搭

リには直接アクセスすることができず，メインメモリとビ載したアーキテクチャとなっている．中国の Tianhe-2 はデオメモリ間でPCI Expressバス（Gen2の帯域は8 GB/sec,

Xeon Phi を約５万個搭載するスパコン（ピ ー ク 性能 55

Gen3 では 16 GB/sec）を介したデータ通信を行うことがで

PFlops）で，Top500 ランキングにおいて３期連続（１年きる．GPU スパコンの各ノードには，２個程度の CPU と半）で世界一になっている［３］．

GPU を搭載したビデオカードが１〜４枚程度搭載されて

１．
３ スパコンの電力効率いることが多く，さまざまなメモリが階層的構造になっていて，それぞれデータ転送の帯域が大きく異なる点が特徴

世界トップクラスのスパコンが PFlops を超えてきたあたりからエクサスケール（Exa-Flops）のスパコン開発の議的である．

論が始まった．当時のスパコンを１０００倍にスケールするGPU スパコンでは，千〜万ノードがさまざまなトポロと，スパコンの消費電力は世界最大級の発電所の電力に匹

ジーの高速なインターコネクションで接続されている．

敵してしまい，スパコンにおいても低消費電力化が最重要

TSUBAME2.5では１ノードに２個のIntel Xeon X5670と３

課題と認識されるようになった．理化学研究所の「京コン個のTesla K20Xが搭載され，
１４０８ノードがQDR InfiniBandピュータ」の消費電力は 13 MW，Tianhe-2 の消費電力は約により Fat Tree 型で接続されている．

18 MW であり，今後も消費電力は 20 MW 程度が限界だと考えられていて，エクサスケールのスパコンを完成させる

２．GPU コンピューティングには，今後数年間の技術革新でプロセッサの消費電力だけ

２．
１ 軽量・超多スレッドの並列計算
前節で述べた通り，１つのGPUには２６８８個もの演算コア

でなく，冷却を含めたスパコン・システムの消費電力を２０が搭載されており，それらを効率よく利用することにより分の１に下げなければならない．
モバイル IT 機器では搭載する CPU の低消費電力化が進

GPU の高い計算性能を引き出すことができる．即ち，並列

められ，ダークシリコンと呼ばれる消費電力の制約からシ計算が必要になる．SMX の中の各 CUDA コアはスレッド

リコンチップ上で電力を供給してオンにできないエリアのという単位で計算を実行するが，詳細には３２個の CUDA

問題解決に GPU の利用が進められている．GPU は消費電コアが同一の演算命令を実行する必要があり，SIMD（Sin-

力当たりの演算性能が高く，スパコンにおいてもアクセラ

gle Instruction Multiple Data）計算となる．例えば，メモリレータとしての導入が進んでいる．もう一つのスパコンの上のデータの数値は違ってもよいが，３２個の CUDA コアがランキングに Green500［４］があり，Top500 のランキング指一斉にメモリアクセスを行い，レジスタにデータを保存標である Linpack スコアを消費電力で割った Flops/W を指し，同一の演算を行う処理になる．SMX に対しては２０４８個標としている．２０１４年６月のGreen 500ランキングは，東京までのスレッドを同時に実行命令することができ，投入さ工業大学・学術国際情報センターの油浸冷却を取り入れた

れたスレッドをスケジューラが３２個単位に分割し，メモリ
７５６

Commentary

Large-Scale Numerical Simulations Based on GPU Computing

T. Aoki

ストコードには GPU のビデオメモリ上への配列確保やメ
インメモリとビデオメモリ間のデータ転送などの API と
GPU の実行内容を記述したカーネル関数のコール等が記述される．デバイスコードにはGPUカーネル関数の内容が書かれ，その内容がスレッドとして各 CUDA コアで実行される．
CPU を使ってメインメモリ上の配列 A と B に対して，n個の要素を B から A にコピーする C 言語のプログラム１
図２

NVIDIA Kepler コア GPU のアーキテクチャ．

は，アクセスと演算処理を同時に実行させるなどの効率的な処

double *A, *B;

理が行われる．GPU 全体に対しては，物理的に存在する

int

２６８８個のコア数をはるかに上回る最大で１０２４×１０２４×６４×i;

for(i = 0; i< n; i++) {１０２４個のスレッドを投入することができる．１つのスレッ

A[I] = B[i];ドの処理内容は１回のメモリアクセスと四則演算のような}軽い処理でもよく，最低でも数万スレッド以上ないと高い実行性能は期待できない．スレッドの内容は条件分岐など

プログラム１を含んでもよく，１つのGPUに対しては同一のプログラム

CPU の逐次コード．

が動く SPMD（Single Program Multiple Data）となっていと簡単に書ける．一方，GPU のビデオメモリ上に確保した

る．

配列 A と B に対して，n 個の要素を B から A にコピーする

複数 GPU を使う場合は，ノード内ならば Open MP など

GPU のカーネル関数に対する CUDA のプログラムは，を利用して複数 GPU を制御し，複数ノードになる場合はMPI ライブラリを用いてそれらが階層的な SPMD として__global__ void array_copy(double *A, double *B)

実行される．

{

２．
２ GPU を使うためのプログラミング

int i = blockDim.x*blockIdx.x + threadIdx.x;

２．
２．
１ CUDA プログラミング

B[i] = A[i];第１章でも少し述べたが，GPU を計算に用いるために}は，殆どの場合は CUDA でプログラミングする必要があ
る．CUDA は C 言語または FORTRAN 言語に GPU コン

プログラム２

ピューティングのための言語拡張や API（Application Pro-GPU のカーネル関数コード．

gramming Interface）やコンパイラなどを含めた統合開発となる．CUDA の C プログラムでも，CPU のプログラム１

環境である．FORTRAN コンパイラのみ有償であるが，そと同じようにポインタや配列を使うことができる．プログ

れ以外は NVIDIA から全て無償で提供され，Linux，Win-

ラム２の内容がスレッドとして１つの CUDA コアで実行

dows，Mac などのほとんどのプラットホームに対して準される．ホストコードに次のようにカーネル関数の実行指備されている．また，サンプル・コードや解説ドキュメン

令が記述される．

トも多く，初心者に対しても敷居が低い［１］．

array_copy<<<n/256, 256>>>(A, B)

CUDA のプログラミングについて概略を図３に説明す（１）

る．全体では CPU での実行を記述するホストコードと“<<<a，b>>>”の中の２つの引数がGPUに投入するスレッ

GPU の実行を記述するデバイスコードから構成される．ホド数を階層的に指定している．第２引数は b 個のスレッドが１つのブロックを構成することを指定し，第１引数の aはブロック数を指定している．したがって，a×b 個のス
レッドが GPU で実行されることになる．a，b が整数の場合
は１次元的にスレッド を 管 理 す る こ と に な る が，２次元，３次元ベクトルとして巨大な総数のスレッドを指定することができる．ブロック内のスレッドは同一 SMX の中で処理されるが，どのブロックがどの SMX で処理される
か，ブロック中のどのスレッドがどの CUDA コアで実行されるかは指定できない．
（１）の例では，カーネル関数コールにより１ブロック内に２５６個のスレッドを配置し，n/256
個のブロックが投入されるので総計 n/256×256＝n 個のス

図３

レッドが実行される．カーネル関数内ではスレッド毎に違GPUで計算を実行するためのホストコードとデバイスコー
ドの関係．

う値が入るビルトイン変数を宣言なしで使うことができ
７５７

Journal of Plasma and Fusion Research Vol.90, No.12 December 2014

る．プログラム２の中で使われている blockDim.x にはブロック中のスレッドの数が入っており，（１）でカーネル関double *A, *B;

数コールを行った場合には blockDim.x＝２５６となる．ビint

ルトイン変数 blockIdx.x にはそのスレッドが何番目のi;#pragma acc kernels!

ブロックの中で実行されるかの値が入り，0〜n/256−1 の#pragma acc loop independent!

どれかの整数になる．threadIdx.x にはそのスレッドがブロック内の何番目であるかが入り，０〜２５５の値になる．

for(i = 0; i< n; i++) {これらのビルトイン変数を利用し，n 個の配列要素に対しA[I] = B[i];て，図４のような対応関係を規定することにより，投入さ}れ る n 個 の ス レ ッ ド に 対 し て，i = blockDim.x*
blockIdx.x + threadIdx.x は 0〜n−1 の全て異なる値プログラム３

OpenACC 指示行の挿入．

を取らせることができる．したがって，プログラム２のカーネル関数を（１）でコールすることにより，n 個のスレッ

定された領域がアクセラレータで実行されるカーネルへ変ドが全て異なる配列要素をコピーすることができ，プログ

換され，GPU で実行される．プログラム３では示していな

ラム１と同じ内容をビデオメモリに対して実行することがいが，data ディレクティブを挿入すると，GPU のビデオメできる．プログラム１はループにより逐次的に同じコピー

モリ上への配列確保，メインメモリ−ビデオメモリ間のを n 回繰り返しているが，プログラム２では，ある１要素

データ転送が制御される．これらはコンパイラ側での処理をコピーするだけというスレッドを n 個投入していて，であり，これまでにさまざまな不具合が解消されてきていループがなくなっている点が興味深い．スレッドの総数はて，現 時 点 で PGI，CAPS，Cray か ら バ ー ジ ョ ン２．
０のCUDA コアより圧倒的に多いので１コアは複数のスレッ

OpenACC 対応のコンパイラが有償で提供されている．

ドを実行するが，その実行順はスレッド・スケジューラに

GPU コ ン ピ ュ ー テ ィ ン グ と い え ば，現 在 の と こ ろ

より決められ，プログラム側で指定することはできない．

NVIDIA の GPU に対してピーク性能に近い性能まで引き

CUDA は GPU の性能向上に合わせ，機能向上を含んだ出すことを目的とする CUDA と，既存のプログラムに対し頻繁なバージョン・アップが行われ，現時点でバージョンて移植のコストを抑えつつ，ある程度の高速化を期待する

６．
０が正式にリリースされている．

OpenACC の二つが中心といえる．

NVIDIA の GPU だけでなく，AMD 社や Intel 社の GPU
や IBM 社の Cell，DSP などに対して，同じプログラムが動３．GPU スパコンによる大規模シミュレーション

作するような標準化（クロスプラットフォーム）として

３．
１ 複数 GPU を用いる計算の問題点

Open CL（Open Computing Language）の仕様が２００９年に単一のグラフィクス・カードに搭載されるビデオメモリ

策定され，KhronosGroup により現在までにバージョン２．
０のサイズはせいぜい 6 GB（最新の Tesla K40 では 12 GB）でが リ リ ー ス さ れ て い る．初 期 は NVIDIA，AMD，あり，より大規模計算を行うには複数のGPUを使った計算Apple，Intel 等の多くの企業が OpenCL コンソーシアムにが必要となる．格子系アプリケーションでは領域分割法が参加し大きな期待がもたれたが，NVIDIA 社が自社の GPU用いられ，分割された各領域を１個の GPU が計算する．分に対 す る OpenCL へ の サ ポ ー ト を 更 新 し な く な ったた

割された領域の境界近傍格子での計算は隣接領域の格子点

め，開発の勢いは鈍化している．

にアクセスする必要があり，袖領域のデータ通信が発生す２．
２．
２ 指示行挿入による GPU の利用（OpenACC）

る．ビデオメモリ上にあるデータのGPU間通信は帯域の狭２０１２年に NVIDIA 社，PGI 社，CAPS 社，Cray 社が協定

い PCI Express バスを介して通信しなければならずメモリ

し，CPU で実行する通常のソースコードに指示行を挿入すアクセスと比較するとかなり遅い．さらにCPU側のメモリるだけでGPUの機能を使えるようにするOpenACCの開発を介して転送する必要があるため，大規模計算においてが始まった．ホスト（CPU）用コード，さまざまなアクセ

GPU 間のデータ通信は大きなオーバーヘッドになる．最新

ラレータ用コードを単一コードとして記述でき，メンテナの NVIDIA の GPU Direct RDMA を利用することで異なるンスが容易になるなどの利点がある．並列化したいループ

ノードのGPUとの通信の手続きが簡単になり，通信速度もに kernels，loop ディレクティブを追加することにより，指２倍近く早くなるが，依然として大きなオーバーヘッドがあることは変わりがない．高い実行性能を得るためには，分割された各領域において境界格子を先に計算し，その終了とともにGPU間データ通信を開始する．それと同時に内側の格子の計算を並行して実行する「計算と通信のオー
バーラップ（図５）」により，GPU 間データ通信の時間を可能な限り隠蔽することが必要である．
３．
２ 500 m 水平解像度を用いた次世代気象計算

図４

気象計算はスパコンを利用する大規模計算の代表的アプ

ビルトイン変数と配列要素の対応．

リケーションである．次期気象予報モデルとして気象庁が７５８

Commentary

Large-Scale Numerical Simulations Based on GPU Computing

T. Aoki

TSUBAME2.0 の全ノードを用いた計算ではピーク性能の１５％となる 600 TFlops,TSUBAME2.5 では 1.14 PFlops の実
行 性 能 が 得 ら れ た．１０，
０８０×１０，２４０×５１２格 子 に 対 し て
４，
０３２個の GPU を用い，新宿や皇居を含む 10 km 四方のエリアを 1 m 格子間隔で計算した．図７は風速分布を多数の粒子を用いて可視化している．高層ビル背後の発達した渦によるビル風や幹線道路に沿って流れる「風の道」，台風の際の被害などが飛躍的な精度で予測できるようになる．
図５

さらに，排ガス，事故やテロによる有毒ガスなどの汚染物

計算と通信のオーバーラップ．

質の拡散も詳細に予測可能となる．
３．
４ 気液二相流シミュレーション

開発している ASUCA に対し，物理過程と力学過程の全てのモジュールのGPU化を行った［５］．力学過程は流体計算

気体と液体が激しく混じり合うような流れ（気液二相そのものであり，双曲型の圧縮性流体方程式を陽的に時間流）は数値流体力学の中でもチャレンジングな研究テーマ

積分しているため，GPU コンピューティングに適していとして知られている．計算時間もかかるため，予てからる．しかし，単純な CUDA による GPU コンピューティン

GPU 計算の導入が望まれていた．GPU コンピューティン

グへのポーティングではなく，さまざまな高速化アルゴリ

グが当初重力多体問題での粒子計算で成功したことが影響

ズムおよび前節で述べた計算と通信のオーバーラップも導し，自 由 表 面 を 含 む 流 れ の 計 算 に 対 し て も，SPH入した．図６は，気象庁で現在の気象予報で実際に使われ（Smoothed Particle Hydrodynamics）法などの粒子法が使ているメソモデル（MSM）のデータを初期条件・境界条件

われてきた．しかし，従来の粒子法はランダムなメモリアと し，TSUBAME 2.0 の 437 GPU を 用 い て GPU 版 のクセス，演算量，計算精度の三つの観点から効率が悪い．

ASUCA により４
７９２×４６９６×４８（水平方向に 500 m×500 m

TSUBAME2.0 を用いて，VOF（Volume of Fluid）法をベーの格子）の格子で，初期時刻２００
９年１０月６日 15 UTC からスとした格子法による気液二相流の計算をマルチ GPU で計算した９時間後の雲の分布を可視化している．気象庁が実行し，高い実行性能を達成することができた．非圧縮性

現在行っている気象予報は水平 5 km および 2 km 格子を用

Navior-Stokes 方程式を半陰解法により時間積分し，移流いているが，その 1/10〜1/4 の格子サイズでの計算が既に項には５次 WENO スキームを使い高精度計算を行った．

TSUBAME 2.0 によりほぼ実時間で可能であることが示さ

気体と液体の界面を記述するためにレベルセット法とカッれ，世界に先駆けGPUスパコンの実用的な運用アプリケー

プルした VOF 法を導入し，毎ステップ再初期化を行ってションに対する有用性が実証された．最近では，計算機いる．

アーキテクチャの急速な進歩に対するソフトウェアの保守性を確保するために，DSL（Domain Specific Language）やフレームワークを導入する研究が盛んに行われている［６］．
３．
３ 1 m 格子による東京都心部の 10 km 四方の気流計算
都市は高層ビルが立ち並び複雑な構造をしており，詳細な気流を解析するためには高解像度格子による大規模気流
シミュレーションが必要となる．数値計算手法は単純なアルゴリズムで大規模計算に適した D3Q19 モデルの格子ボ
ルツマン法［７］を用いた．都市の気流はレイノルズ数が１００
万を超えるような乱流状態になるため，ラージエディ・シ
ミュレーション（LES）［８］のモデルを導入する必要があ
る．現在よく使われている動的スマゴリンスキー・モデル

図６

水平解像度 500 m 格子を用いて次世代メソスケール気象モ
デル ASUCA で計算した雲分布．

図７

1 m 格子を用いて 10 km 四方の範囲の都市部の気流計算（粒子は気流の可視化のために使用）．

［９］では，モデル定数を決定するために各格子点で広領域の平均操作が必要になり，大規模計算には極めて不向きである．本研究では，モデル定数を局所的に決定できるコ
ヒーレント構造スマゴリンスキー・モデル［１０］を格子ボルツマン法に導入することに成功し，大規模な気流の LES
計算を初めて可能にした．実際の建物データに基づき計算対象のエリアを領域分割し，TSUBAME2.0/2.5 の GPU を用いて計算を行った．CUDA を用いてコードを実装し，ここでも３．
１節の計算と通信のオーバーラップを導入し，実行性能をオーバーラップしない場合と比較して３０％以上向上させることができた［１１］．
格子ボルツマン法はメモリ律速の計算手法であるが，
７５９

Journal of Plasma and Fusion Research Vol.90, No.12 December 2014

一番計算時間のかかる Poisson 方程式の解法は，気体とる．導出される方程式は時間空間の偏微分方程式になって液体の密度比が非ゼロ要素に含まれ，みずほ情報総研と共いて，有限差分法や有限要素法などの格子法で解かれるこ同研究している Bi-CGSTAB 法［１２］にマルチグリッド前処とが多い．フェーズフィールド法は複雑な非線形項を多く

理（図８）を導入したライブラリを用いている．!はラプ含むための１格子点あたりの演算量が多く，安定性の観点

ラシアンの離散オペレータ， #は従属変数（ここでは圧から時間ステップを小さく取らなければならない．CPU力），"は右辺のソース項（ここでは速度の発散）を示し，で計算すると時間がかかり過ぎるため，これまでは主に２上付き添え字は格子の粗度を表している．粗度が０は解く

次元計算による解析が行われてきた．大規模並列計算でべき格子解像度であり，１増える毎に格子間隔は２倍，格は，３次元領域分割を行うことにより領域間のデータ通信

子点数は 1/2 になる．各段で残差をソース項として修正ベ量を最小にすることができるが，GPU 計算の場合は $方向

クトルを求めるためのスムーザとしての不完全 ILU 分解をで分割した場合の $!%断面のメモリアドレスが不連続とな行っている．GPU では格子間で順序依存性がある計算はでり，メモリアクセス性能の低下のオーバーヘッドの方が通きないため，ここでは Red and Black 法を用いている．ま

信量の低減より大きくなる．ここでは，%方向，&方向に分た，細かい格子から粗い格子への制限補間（Restriction）

割する２次元の領域分割法を用い，それぞれの分割領域のと粗い格子から細かい格子への延長補間（Prolongation）

計算を各 GPU へ割り当てる．隣接する GPU 間での境界領は単純な平均化処理を行っている．CPU での計算と同じよ域のデータ交換はホスト側のメモリを経由し，次の３段階

うに図８の V サイクルを GPU でも実装し効率的な収束計で構成される．
（１）CUDA ランタイムライブラリによる算を行っている．

GPU から CPU への転送，（２）MPIライブラリによるノード

計算では表面張力や壁との接触角も考慮し，単一気泡の間のデータ転送，（３）CUDA ランタイムライブラリによる上昇を実験と詳細に比較し十分な一致が得られ，図９に示

CPU から GPU への転送を行う．ここでは隣接ノード間です濡れた床へ侵入するダムブレーム問題では，早期に砕波のデータ通信と計算に対し，以下の３つの方法を実装し比が起こり水と空気が激しく混じり合い，壁への衝撃圧が低較を行った．

下する様子が再現されている．九州大学（応用力学研究所（a）GPU-only method：比較のために計算と通信のオー

胡長洪教授）との共同研究において，壁での衝撃圧の時

バーラップを行わず，上記の通信の３ステップを順に間履歴が計算と実験で定量的によく一致することを確認し行う．
（b）Hybrid-YZ method：１つの GPU が担当するサブ領域

た．
３．
５ フェーズフィールド法による合金の樹枝状凝固シを，%方向の境界領域，&方向の境界領域，残りの中心

ミュレーション

領域に分割する．%，&方向の境界領域を CPU で計算

金属材料の組織構造は凝固過程のミクロな結晶構造で決し，中心領域を GPU で計算する．先行研究［１４］では境定される．一方，現実に機械的強度を判定するには，数界も中心領域もGPUで計算したが，ここでは境界領域

mm のマクロなスケールでの解析が必要となる．ミクロなを CPU で計算し GPU 関数を分割することなく通信を凝固ダイナミクスを解明するために，非平衡統計力学から隠蔽する．ただし，全計算領域サイズによって CPU導出されるフェーズフィールド法［１３］が近年注目されていによる通信と境界領域の計算にかかる時間が中心領域の GPU 計算時間よりも長いことがあり，CPU がボト
ルネックとなることが考えられる．
（c）Hybrid-Y method：（b）と同じように CPU 計算も利用するが，&方向の境界領域（$!%断面）はカーネル関数
を分割して GPU で計算する．境界領域を担当する
GPU カーネル関数は，グローバルメモリの連続アクセスであるため計算効率がよく，また通信バッファへ
データを詰め替える必要もない（図１０）．

図８

図９マルチグリッド前処理の V サイクル．

VOF 法ベースの格子法による気液二相流シミュレーション．

図１
０ Hybrid-Y method のダイアグラム．

７６０

Commentary

Large-Scale Numerical Simulations Based on GPU Computing

T. Aoki

TSUBAME の 各 ノ ー ド に は 3 GPU と 2 CPU sockets

用と違い，砂や粉などの粉体現象は粒子間の接触による反（12CPU コア）が搭載されている．このため CPU による境発 や 摩 擦 力 の モ デ ル（離 散 要 素 法：Discrete Element

界 領 域 の 計 算 に 1 GPU あ た り 4CPU コ ア を 割 り 当 て，

Method（DEM））でシミュレーションされる．実現象と同OpenMP を用いた並列計算を行っている．TSUBAME2.0

じ程度の粒子サイズで計算することにより，疎視化モデルの 複 数 GPU
（M2050）を 用 い（a）GPU-only method，（b）では表現できない現象を計算することが可能になるが，計Hybrid-YZ method，（c）Hybrid-Y method の３つの実装の算規模が膨大になり単一 GPU ではメモリが足りず複数実行性能について述べる．フェーズフィールド法により Al

GPU による計算が必要になる．

-Si 合金の凝固成長の計算を行った．GPU コードの実行性

重力多体問題などと異なり DEM 計算はメモリ律速であ能の評価において，GPU コードの浮動小数点演算数を直接り，粒子の位置や速度などの従属変数はGPUのビデオメモ

測定することは困難である．ここでは全く同じ結果を出す

リ上に乗せておく必要がある．また，接触による相互作用

CPU コードに対して PAPI（Performance API）を用いて実であるため，粒子分布を空間領域で分割し，分割された領測した浮動小数点演算数を元に GPU の実行時間から実行

域内の粒子数を一定にすることで GPU 間の計算負荷を均性能を評価している．

一にして並列計算の効率を上げることができる．しかし，

図１１の強スケーリングは，問題を固定してGPU数を増や粒子の空間分布は時間とともに大きく変化し，初期に均一

して行くとき の 実 行 性 能 の 向 上 を 示 し て い て，前 述のな負荷であったとしても時間とともに負荷バランスは大き（a），（b），（c）の３つの実装の違いによる特徴が表れていく崩れる．そこで，図１４のように領域境界を移動させるス

３

３

３

１節の
る．全体の計算領域を５
１２ ，１０２４ ，２０４８ とした．３．

ライスグリッド法を導入することにより，常に領域内の粒ように計算と通信のオーバーラップを導入した（b）Hybrid-

子数をほぼ一定に保つ動的負荷分散を行うことができる

YZ method，（c）Hybrid-Y method では GPU 数が少ないと［１６］．

きには通信を隠蔽できている．GPU 数 の 増 加 と と もに

CAD データで表現される物体形状を符号付距離関数のGPU の計算時間が短くなるため，ある段階で計算時間より通信時間の方が長くなり，もはや通信を隠ぺいすることができなくなる．
（c）Hybrid-Y method では，GPU 数を増やした時に CPU の計算時間が隠ぺいの足を引っ張る部分が大幅に改善され，広い範囲のGPU数で最適化を導入していない（a）GPU-only method と比較して実行性能が大幅に改善されている．
弱スケーリング（図１２）は，1GPU あたりで実行する問題
規模を一定にして GPU 数を増やしたときの性能評価である．この測定では，1GPU あたりの計算格子サイズを単精度計算では４０９６×１２８×１２８，倍精度計算では４０９６×１２８×６４とした．弱スケーリングの結果が得られ，（c）Hybrid-Y

図１
２ 樹枝状凝固シミュレーションの弱スケーリング．

method では，4000 GPU を利用して４０９６×６５００×１０４００計
算 格 子 に 対 し て 行 っ た 計 算 で，2.0 PFlops（GPU:1.975
PFlops，CPU:0.025 PFlops）という極めて高い実行性能を達成した．この成果により２０１１年に ACM ゴードンベル賞を受賞している［１５］．図１３は，Al-Si 合金の凝固成長シミュ
レーションのスナップショットを示している．
３．
６ 動的負荷分散を導入した粉体シミュレーション
重力やクーロン力のような長距離力による粒子間相互作

図１
３ Ai-Si 合金の樹枝状凝固過程のシミュレーション．

図１
４ スライスグリッド法による粒子計算の動的負荷分散．

図１
１ 樹枝状凝固シミュレーションの実行性能（強スケーリング）．

７６１

Journal of Plasma and Fusion Research Vol.90, No.12 December 2014

等値面で近似することによりポリゴンのエッジの特異点を排除し，物体との接触判定の計算負荷を大幅に低減することができる．
TSUBAME2.5 の 64GPU を用いて粉体の撹拌シミュレー
ションを行った．図１５は５００万個の粒子を用いているが，256 GPU により１億３０００万個の粒子計算に対しても効率よ
く計算が行われている．

４．Plasma シ ミ ュ レ ー シ ョ ン へ の GPU コ ン
ピューティングの適用
４．
１ MHD シミュレーション
前節で示した気象計算の力学過程は圧縮性流体計算を行うために空間を格子で離散して Navior-Stokes 方程式を解くもので，MHD 方程式の離散化とほぼ同じである．MHD
図１
６ 太陽風と惑星の相互作用により発生する衝撃波．

シ ミ ュ レ ー シ ョ ン の 陽 的 時 間 積 分 は 容 易 に GPU コ ン
ピューティングを行うことができる．Maxwell 方程式の電磁波モードについては陽的時間積分であるので問題ないする必要がある．また，ベクトル型計算機でも工夫が必要

が，静電場をポアソン方程式で解く場合は楕円型方程式にになった粒子分布から格子点上の電荷密度や電流密度を求なり，大規模疎行列計算が必要になる．３．
４節で説明した半

める gathering 処理は足し合わせがどうしても逐次処理に陰解法のための圧力のポアソン方程式は静電場の計算と全なり，並列処理を行う有効なアルゴリズムがない．CUDA

く同じであり，拡散モデルによる輻射輸送は放物型方程式には排他的に他のスレッドの書き込みを止め，逐次処理をであり楕円型方程式の計算に準じ，疎行列に対するクリロ

可能にする ATOMIC 処理が可能であるが，通常は非常にフ部分空間法の反復計算や有効な前 処 理 を 同 じ ように時間のかかる処理となり高速計算には不向きである．しかGPU 計算で使うことができる．Vlasov-Maxwell 方程式やし，荷電粒子の gathering は Cell 内の局所的な足し合わせFokker-Planck 方程式を解く Plasma Kinetics シミュレー

であり，CUDA のブロック内で完結させることができるた

ションについても，ほとんどの部分が双曲型方程式であめ，足 し 合 わ せ を SMX 内 の 高 速 な 共 有 メ モ リ 上 でり，GPU コンピューティングを行う上での問題はない．電ATOMIC 処理を用いて行うことにより高速に計算を行う

離や励起などの原子過程は化学反応とほとんど同じ密行列ことができる［１９，２０］．

計算になり GPU の得意とする計算になり，CUDA の BLAS

５．国内の GPU スパコンの利用についてや Lapack ライブラリが利用できる．
筑波大学のグループは日本原子力研究開発機構のプラズ

GPU コンピューティングをプラズマ・核融合分野のシ

マ乱流 を 解 析 す る た め の５次 元 Vlasov コ ー ド GT5D のミュレーションに適用することで，これまでできなかった

GPU 化を行っている［１７］．また，太陽風と惑星の相互作

規模の計算を短時間で実行し，新しい現象の発見や理論の用に対する大規模な MHD シミュレーション（図１６）も行

解明につなげる可能性がある．CUDA のプログラミングやわれている［１８］．

OpenACC は確かに取っ付き難いところがあるが，その後

４．
２ PIC シミュレーションはそれほど難しいものではない．ほんの少しの努力で最初の障壁を乗り越えれば，得るものは大きいと信じている．

PIC（Particle-in-Cell）シミュレーションも粒子計算に
GPU を適用することは容易である．格子点上で電磁場を計

国 内 に は 東 京 工 業 大 学・学 術 国 際 情 報 セ ン タ ー の算するのは負荷が低く，大規模計算を行うには３．
６節のよ

TSUBAME2.5 以外にも，京都大学学術情報メディアセン

うな何らかの粒子の空間分布に対する動的負荷分散を導入

ター，九州大学情報基盤研究開発センターなどにも小規模の GPU を搭載したスパコンがある．学術的課題に対しては，それらを学際大規模情報基盤共同利用共同研究拠点の公募型共同研究課題や「京コンピュータ」以外の HPCI（革新的ハイパフォーマンス・コンピューティング・インフラ）構成機関が提供する計算機資源を使う課題として採択されれば一定量の計算資源が無償で利用できる．また，産業利用についても，東京工業大学・学術国際情報センター
では無償のトライアルユース，有償の共同利用サービス（成果公開／非公開）を行っており，GPU を利用できる機会は確実に広がっている．

図１
５ 離散要素法による粉体の撹拌シミュレーション．

７６２

Commentary

謝Large-Scale Numerical Simulations Based on GPU Computing

IEEE Computer Society, New Orleans, USA (2014).
［７］X. Wang and T. Aoki, Parallel Computing, 37, 521 (2011).
［８］M. Lesieur et al., Large-eddy simulations of turbulence (Cambridge University Press, New York, 2005).
［９］M. Germano et al., Phys. Fluids A: Fluid Dynamics, 3, 1760
(1991).
［１０］H. Kobayashi, Phys. Fluids 17, 045104 (2005).
［１１］小野寺直幸他：情報処理学会ハイパフォーマンスコン
ピューティング研究会主催HPCSシンポジウム（２
０
１
３）
．
［１
２］H.A. van der Vorst, SIAM J. Sci. and Stat. Comput. 13, 631(1992).
［１３］R. Kobayashi, Physica D, Nonlinear Phenomena, 633-4,
410 (1993).
［１４］小川 慧他：情報処理学会論文誌コンピューティングシ
ステム 3, 67 (2010).
［１５］T. Shimokawabe et al., in Proceedings of the 2011 ACM/IEEE International Conference for High Performance
Computing, Networking, Storage and Analysis, SC'11,
IEEE Computer Society, Seattle, USA (2011).
［１６］都築怜理他：情報処理学会論文誌コンピューティング
システム 9, 82 (2013).
［１７］N. Fujita et al., Proc. of Int. Workshop on Parallel and Distributed Scientific and Engineering Computing 2014
(PDSEC2014) (with IPDPS2014), Phoenix, USA 2014.
８］Un-Hong Wong et al., Comput. Phys. Commun. 185, 1901
［１(2014).
［１９］X. Kong et al., J. Comput. Phys. 230, 1676 (2011).
［２０］V.K. Decyk and T.V. Singh, Comput. Phys. Commun. 182,
641 (2011).

辞

本解説で紹介した計算は TSUBAME2.0/2.5 で実施したもので，東京工業大学・学術国際情報センターに深く感謝の意を表する．本研究の一部は科学研究費補助金・基盤研
究（S）課 題 番 号 ２６２２０００２「も の づ く り HPC ア プ リ ケ ーションのエクサスケールへの進化」および基盤研究（B）課
題番号 ２３３６００４６「GPU スパコンによる気液二相流と物体の相互作用の超大規模シミュレーション」，科学技術振興機構 CREST「ポストペタスケール高性能計算に資するシ
ステムソフトウェア技術の創出」，学際大規模情報基盤共
同利用・共同研究拠点，および革新的ハイパフォーマン
ス・コンピューティング・インフラから支援をいただいた．記して謝意を表す．

参考文献［１］青木尊之，額田 彰：はじめての CUDA プログラミング（工学社，２００９）．
［２］NVIDIA: CUDA Programming Guide 6.0 (2014).
［３］http://www.top500.org/［４］http://www.green500.org/［５］T. Shimokawabe et al., in Proceedings of the 2011 ACM/IEEE International Conference for High Performance
Computing, Networking, Storage and Analysis, SC'11,
IEEE Computer Society, New Orleans , USA (2010).
［６］T. Shimokawabe et al., in Proceedings of the 2014 ACM/IEEE International Conference for High Performance
Computing, Networking, Storage and Analysis, SC'14,!!!!!!!!!!!!!!!!!

あお

き

たか

ゆき

!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!

青 木 尊 之東京工業大学学術国際情報センター教授．
助手の時代まで，慣性核融合の理論・シミュレーションの研究に従事．最近は大規模数値流体力学，計算力学，ハイパフォー
マンス・コンピューティングを専門とし，登録者数１，０００名を超えるGPUコンピューティング研究会を主宰し，これまでに１７回のハンズオン講習会を開催する等，CUDA プログラミ
ングの普及も行っている．２
０
１
１年に ACM コードンベル賞を受賞している．

T. Aoki!!!!!!!!!!!!!!!!!!!!!!!!!!

７６３