― 54 ―GPU を用いた LGA 法の高速化
高橋直也＊・竹下大樹Acceleration of LGA method using GPU
Naoya TAKAHASHI ＊ and Daiki TAKESHITA（平成23年11月25日受理）
LGA method supposes a fluid to be the set of the particle, and it expresses a flow of fluid
from interaction between particles. The interaction is calculated in the whole system all at once.
Therefore, LGA method is theoretically suitable for parallel computation. Then, it attemptedby using GPU that the parallel processing ability is high for the speed-up of the fluid simulation
program by the LGA method. CUDA was used for the programming of GPU. As a result,as compared to programs using only CPU, the processing speed was improved 6-7 times.

1.

緒論

格子ガスオートマトン法（以下 LGA : lattice gascellular automaton）法は流体をシミュレーションする手法の一つである。この手法では平面を格子状に分割し，粒子の流れを格子線上に制限する。そして，粒子間で質量と運動量を保存する衝突を繰り返し，流体の速度場を解析する。この相互作用は系全体で一斉に計算されるため，原理的に並列処理に適している。また，近年，画像処理を専門とする補助演算装置である GPU を画像処理以外の目的に応用する GPGPU という技術が注目されている。GPU は一度に大量のデータを処理することに特化しており，CPU よりも多くの演算ユニットを搭載している。本研究ではこの GPU を用いて LGA 法による流体シミュレーションの高速化を図り，CPU のみの場合と比較，検討をすることを目的とする。なお，
GPU を用いたプログラムの作成には，NVIDIA 社の提供する GPGPU 専用統合開発環境である CUDA（Compute Uniﬁed Device Architecture）を用いた。
2.

れる。また，各粒子には時間と境界について初期条件が与えられており，各粒子は時間ステップごとに格子線上を通り，最近接の格子点まで単位速度で移動する。格子点上では，別の粒子と衝突・散乱して速度方向が変わる。そして，この移動・衝突・散乱のプロセスをあらかじめ決めた衝突側に従ってすべての格子点で一斉に行い，平衡状態になるまで繰り返す。

図1

三角格子状に分割された平面

粒子の移動方向のパターンは図 2 のように格子点上の粒子を中心として計 7 種類で表され，LGA 法によるシミュレーションはこの 7 種類の状態の粒子の有無により行われる。

LGA 法

LGA 法では，平面は図 1 に示すように正三角形を敷き詰めた格子に分割され，単位質量の粒子は格子点上に存在しており，その移動は格子線上に限ら＊秋田高専専攻科学生

図2

一つの格子点上での粒子の移動方向
平成24年２月

― 55 ―GPU を用いた LGA 法の高速化

本研究では衝突側に FHP モデルを用いた。FHP
モデルにおいて，考慮される衝突パターンを表 1 に示す。なお，白丸は格子点上に停止している粒子である。
表1

FHP モデルの衝突例

図3

CUDA でのプログラムの流れを集めたものをグリッドと呼ぶ。プログラム中ではこのブロックとグリッドという単位も加え，スレッドを管理していく。なお，ブロックは最大 3 次元の配列として，グリッドは最大 2 次元の配列として定義できる。スレッド，ブロック，グリッドの関係を
図 4 に示す。

3.

CUDA

3.1 プログラムの構成
CUDA によるプログラムは大きくホストプログ
ラムとカーネル関数に分類される。ホストプログラ
ムは，CPU 上で実行されるプログラムであり，C 言語によるプログラムを基本として，GPU を用いるにあたってデータ転送用のメモリ領域の確保やカーネル関数の定義，呼び出しも行う。一方，カーネル関数とは GPU 上で実行されるプログラムである。
このプログラムはホストプログラムから呼び出されることで実行される。CUDA でのプログラムの流れを図 3 に示す。
3.2 スレッドによる管理
CUDA では，プログラムをスレッド単位で管理し，大量のスレッドを同時に実行させる形で GPUによる並列処理を行う。さらに，このスレッドを集めたものをブロックと呼び，同じサイズのブロック
秋田高専研究紀要第47号

図4スレッド，ブロック，グリッドの関係

具体的なスレッドの管理はスレッド固有の ID を用いて行う。カーネル関数を定義するとき，ホストプログラムではグリッドとブロックのサイズを決定する。そして，CUDA ではこのサイズに従って，以下の型の変数が自動的に用意され，そのスレッドに対応した値が格納されている。
（1）threadIdx 型「threadIdx.x」や「threadIdx.y」と記述することで，ブロック内でのスレッドの ID を取得することができる。
（2）blockIdx 型
threadIdx 型と記述の仕方は同じで，グリッド内でのブロックの ID を取得することができる。
（3）blockDim 型

― 56 ―高橋直也・竹下大樹

ブロックの各次元の大きさを取得することができる。
（4）gridDim 型
グリッドの各次元の大きさを取得することができる。
これらの変数を用いることで ID を計算により求めることができる。例として図 5 のグリッドを考える と， 各 ス レ ッ ド の ID は blockIdx.x×blockDim.x＋threadIdx.x により求めることができる。図中の灰色のスレッドの場合，ID を計算するために必要な各変数の値はblockIdx.x：1
blockDim.x：4
threadIdx.x：2となり，スレッドの ID はblockIdx.x ＊ blockDim.x ＋ threadIdx.x＝1＊4＋2＝6となる。

（2）グローバルメモリ
全スレッドから利用可能な大容量のメモリ領域であるが，アクセス速度は低速である。
（3）シェアードメモリ
オンチップメモリであり，レジスタと同様に高速なメモリである。ただし，アクセス可能なのは同一ブロック内のスレッドのみである。
（4）コンスタントメモリ，テクスチャメモリ
デバイス上のキャッシュ機能を持つメモリだが，GPU 側では読み込み専用のメモリとなる。
これらのメモリモデルを図 6 に示す。

図6
図5

一次元スレッドのブロックと一次元グリッド

3.3 CUDA でのメモリ
CUDA において，GPU 側で使用できるメモリは以下の種類がある。
（1）レジスタ，ローカルメモリ
スレッド固有のメモリで，レジスタは GPU チッ
プ内に実装されているオンチップメモリ，ローカルメモリは GPU チップ外のデバイスメモリに配置されているオフチップメモリである。よって，レジスタの方がローカルメモリよりもアクセス速度は高速である。カーネル関数内のローカル変数を保持のために用いられ，レジスタのみではメモリが不足する
場合，ローカルメモリをローカル変数の退避領域として用いる。

CUDA でのメモリモデル

3.4 メモリへのアクセス
プログラミングを行うときメモリへのアクセス方法によって，その処理速度に違いができる。特徴的なのものは以下の二つである。
（1）コアレスアクセス
グローバルメモリへのアクセス方法に関するものである。CUDA では一度に32スレッドが実行され，これをワープと呼ぶ。グローバルメモリにアクセスするときスレッドが連続したメモリに同時にアクセスするようにすると，ハーフワープに相当する16スレッド分のデータを一度の処理で転送することができる。これがコアレスアクセスである。そのためには①

データ型のサイズが32，64または128bit
平成24年２月

― 57 ―GPU を用いた LGA 法の高速化②

ハーフワープの先頭スレッドが64byte または128byte 境界にアライメントされたアドレスにアクセスする。

コンフリクトを回避するようにメモリアクセスすることが必要となる。シェアードメモリへのアクセスの例を図 8 に示す。

といった条件があるが，コアレスアクセスにしない
場合，アクセス命令の分だけ転送処理が必要となる上に通常の演算命令が 4 ～数十サイクルかかるのに対し，グローバルメモリからの転送は400～600サイクルもかかるため，処理速度は大きく低下する。そのため，グローバルメモリを使用するときは，コアレスアクセスにすることが必要となる。グローバルメモリへのアクセスの例を図 7 に示す。

図84.

図7

グローバルメモリ（float 型）へのアクセスの例（2）バンクコンフリクト
シェアードメモリへのアクセス方法に関するものである。CUDA ではシェアードメモリをバンクという単位に16分割する。そして，ハーフワープの16スレッドが別個のバンクにアクセスすることで全てのバンクを使い，並列処理を行う。しかし，複数のスレッドが同じアドレスにアクセスする場合はアクセス先バンクの衝突が起こる。これがバンクコンフ
リクトである。バンクコンフリクトが起こると，一つのバンクを複数のスレッドから読み込むことはできないため，そのスレッドの数だけシェアードメモリへのアクセスを行うことになり，効率は低下する。
よって，シェアードメモリを用いる場合は，バンク
秋田高専研究紀要第47号

シェアードメモリへのアクセスの例

GPU 上での LGA 法の実装

本 研 究 で は NVIDIA 社 の GeForce9800GTX＋ を用いた。デバイスの情報を表 2 に示す。表におけるコアの数とは GPU の演算ユニットの数であり，ストリーミングプロセッサ（SP）の数を表す。SP をグループ化したものがマルチプロセッサ（SM）である。よって，SM 一つあたり 8 個の SP を持つことになる。この SP が CUDA でのスレッドに相当し，同じ SM 内の SP は 4 サイクルにわたって同じ命令を繰り返す。このため CUDA では一度に32スレッドを実行することになっている。
表 2 のデバイスの情報を踏まえて，本研究では22464（縦128×横208）個の格子点をシミュレーションの対象とし，一つの格子点に一つのスレッドを対応させた。この格子点の数は CPU のみで行っていたプログラムのサイズが416×256であり，その比率を残したままグリッドの最大サイズを超えない格子点数を確保した。また，1 ブロックのスレッド数は64（縦 8×横 8）とした。ここで，CPU のみを用いた LGA 法のフローチャートを図 9 に示す。
GPU には図 9 の灰色で示された処理を担当させ― 58 ―
高橋直也・竹下大樹

表2

GeForce9800GTX ＋のデバイス情報

マルチプロセッサの数

16

コアの数

128

ワープ

32

グローバルメモリの合計

511MB

ブロック当たりの最大スレッド数512

ブロックの各次元での最大サイズ

512×512×64

グリッドの各次元での最大サイズ

65535×65535×1

ブロックあたりのシェアードメモリの合計量

16KB

クロック周波数

1.84GHz

図10

CUDA を用いた LGA 法のフローチャート

図11

図9

分割された平面と格子点

CPU のみを用いた LGA 法のフローチャート

る。描画処理を行うまでは CPU には結果を転送せず，GPU とグローバルメモリ間のみで粒子情報の更新を行う。GPU を用いた LGA 法のフローチャー
トを図10に示す。GPU を用いた場合，CPU 側で初期値を設定した後，転送用のメモリ領域を用いてそのデータを GPU へ転送する。その後，カーネルを実行することで各格子点に対応するスレッドがそれぞれ粒子の移動から粒子情報の更新までを行い，結果を CPU へと転送する。CPU 側では描画処理を行うまでに必要な繰り返し回数だけカーネルを実行したかを判断し，再びカーネルを実行するか，グローバルメモリから結果を受け取って粒子の流れを描画する。終了条件は，シミュレーション全体で十分な回数だけ計算を繰り返したかである。
さて，グローバルメモリ上で粒子の情報の更新を行うとしたが，これは GPU 側の粒子の移動処理のためである。説明のため，図11に示す格子上に分割

図12

図11の格子点の並びに対応した配列された平面を考える。番号はその格子点を表すものとする。また，図11の格子点に対応した配列 cell を図12に示す。
粒子の移動処理では，各格子点の周囲 6 点の粒子情報からその格子点に向かってくる粒子の数を判断する。例えば，図10の格子点 7 では図13（a）のように格子点 0，1，6，8，12，13の情報を参照する。また，境界に位置する格子点では平面の格子点の並びが連続するものとして考える。つまり，格子点 0 の場合は図13（b）のようになる。
今，2×2 のブロック 6 つからなるグリッドを考え，各ブロックのシェアードメモリに対応する配列
平成24年２月

― 59 ―GPU を用いた LGA 法の高速化

図15
図13

対象格子点とアクセスする配列の位置関係

参照する格子点

図16
図14

各ブロックのシェアードメモリのデータの一部をコピーする。すると，各ブロックのシェアードメモリには図14のようにデータが格納される。
図14のシェアードメモリのデータを用いて粒子の移動処理を行うとする。しかし，シェアードメモリは同一ブロック内のスレッド以外はアクセスできないという特徴がある。図13の場合，参照するデータが異なるブロックのシェアードメモリに存在するためアクセスすることが出来ず，処理を行うことができない。ブロックのサイズを大きくしたとしても，いずれブロックを越えてアクセスする必要が出てくる。よって，シェアードメモリを用いる方法は粒子の移動処理には適さないことがわかる。
次にグローバルメモリを用いて処理を行うとする。図12のデータ構造のままグローバルメモリに格納し，各スレッドからアクセスする場合，スレッドの ID と配列の添え字が対応しているとすれば，格子点 7～10，13～16のアクセスする配列の位置関係は図15となる。なお，tid は中央の格子点に対応するスレッドの ID である。
しかし，境界上に位置する格子点 0～6，12，13および17～23の場合は図15のように一律に表すことができない。例えば格子点 0 の場合は図16（a），格子点23の場合は図16（b）となる。
格子点ごとにアクセスの記述を変えれば処理は行えるが，その場合は分岐処理となる。GPU では単一ワープ内のスレッドが同じ命令を実行するときに秋田高専研究紀要第47号

一律に表せない位置関係最高効率が得られる。単一ワープ内のスレッドで分岐処理が起きた場合は各分岐を順に実行し，分岐先の異なるスレッドはその間処理を行わない。LGA法では衝突判定のときに分岐処理が必要である。シミュレーションのプロセス上この処理は取り除くことができないので，この処理以外での分岐処理はできるだけ無いこと望ましい。そこで，図17のようにあらかじめ図12のデータの周囲に境界上の格子点が必要とするデータを並べたものを用意し，グローバルメモリに転送することにする。
このデータを用いて図17の周囲を除く配列に粒子の移動処理を適用する。すると，格子点のアクセスする位置関係は図18のように一律に表すことができる。

― 60 ―
高橋直也・竹下大樹

図17

図18

新たにデータを配置した配列

一律に表した新たな配列の位置関係

セスが考えられる。粒子の情報はグローバルメモリ上で更新され，シミュレーションはそのデータを基に繰り返している。通常の演算命令よりも多くのサイクルを必要とするグローバルメモリへのアクセスが多いことで，スレッドでのデータ取得までの待ち時間が発生し，速度の低下に繋がっていると考えられる。データへのアクセスの関係で本研究ではグローバルメモリを用いたが，この結果から高速化のためにはグローバルメモリよりもシェアードメモリを使うことが必要であるといえる。ただし，容量は小さく，アクセスに関する制限もあるためプログラムの並列処理化やデータ構造の改善も必要になるといえる。データの構造に関し て は AoS（Array of Structures） で は な く SoA（Structure of Arrays）を用いることが挙げられる。
AoS とは構造体を要素とする配列のことであり，例を図19に示す。一方，SoA は本来関連性がある構造体のメンバを個々に配列にまとめ，それらを再び配列とするものである。例を図20に示す。なお，図19，20の x，y，z は各座標のデータを表すとする。

位置関係を一律に表せたことにより，各スレッドは同じ命令で処理を行う。そのため，分岐処理無しに粒子の移動処理のためのデータへのアクセスを行えるようになる。よって，本研究では粒子情報の更新はグローバルメモリ上で行うことにした。
5.

結果と考察

図19

AoS の例

図20

SoA の例

CPU のみを用いた場合と GPU を用いた場合のシミュレーションの所要時間についての結果を表 3 に示す。
表3実行モデルCPU と GPU での実効結果
処理時間（秒）

実行速度（倍）CPU のみGPU 使用FHP- Ⅰ200

326.3

FHP- Ⅱ

238

356.8

FHP- Ⅲ261

377.1

いずれの衝突パターンでも実行速度が 6～7 倍になっており，シミュレーションの高速化を実現することが出来たといえる。しかし，CPU がスレッド1 つで計算しているのに対し GPU 側では22464個のスレッドを用いており，理論的には計算ユニット数の多さだけ計算速度は速くなるはずである。この原因としては，まずグローバルメモリへのアクSoA を用いると同じ属性を持つ要素の配列ごとに処理を行えばよく，その処理の内容は同じであるから結果として並列処理を記述しやすくなる。
他の原因としては描画処理時のデータの転送が挙げられる。図10のフローチャートに示すように描画処理を行う場合，一度 GPU 側から CPU 側へと結果を転送する。これがオーバーヘッドとなり，処理速度の低下に繋がったと考えられる。
次に，衝突側で比較すると衝突パターンの多い
FHP― Ⅲが最も速度が速くなっている。これは衝突パターンが多いほど分岐処理の負荷が大きくなり，GPU の方が CPU よりも処理能力において優位にな
平成24年２月

― 61 ―GPU を用いた LGA 法の高速化るためと考えられる。
6.

結言（1）LGA 法に CUDA を適用することで 6 ～ 7 倍のシミュレーションの高速化を実現することができ，GPU の処理能力の高さがわかった。
（2）CUDA で扱うメモリにはコアレスアクセスやバンクコンフリクトといったアクセスに関する制限や容量と速度のトレードオフな関係があり，それぞれの特徴を考えてプログラミングする必要がある。また，扱うデータのデータ型にも注意が必要である。
（3）CUDA ではワープに基づいて処理が行われる

秋田高専研究紀要第47号

ため，ワープを意識してブロック，グリットのサイズの決定する必要がある。
（4）GPU の性能を引き出すにはプログラムの並列処理化や並列処理に適したデータ構造などが必要となる。
参考文献
1） 加藤恭義，光成友孝，築山 洋：「セルオートマトン法―複雑系の自己組織化と超並列処理―」，森北出版株式会社，47P
2） 内田智史：「C 言語によるプログラミング―応用編― 第 2 版」，オーム社