プログラム自動生成とCUDAによる格子ボルツマン法を利用した流体シミュレータの高速化

ネットワークシステム研究室

指導教員坂本　直志

13EC101 後藤　徹哉

近年、PGUを汎用的に使用する事が可能になった事や、計算性能の飛躍的な向上からPGUによる数値流体解析の研究が盛んに行われている。従来並列コンピュータはCPUを使った非常に大掛かりなクラスタで構成されて来た。その最もたる物がスーパーコンピュータである。しかし現在、PGUはその内部に非常に多数のコアを持つ様になり、1つのデバイスそのものが並列コンピュータと言えるまでになっている。こうした状況の中で、従来スーパーコンピュータで計算されていた処理をPGUによって行わせる動きが活発になっている。

本研究では流体の物理シミュレータの高速化の手法として、GPPGUを利用した流体シミュレーションプログラムを動的に生成する事の有用性について評価する。

流体と流体力学

流体とは、水や空気と言った力を加えると変形する連続体の事である。流体は小さな力で形が大きく変化する他体積もまた変化し、流体の粘性や摩擦等もその運動に影響を与える。この流体の運動を研究する為に流体力学と呼ばれる力学の1分野が存在する。

流体力学における流体の運動はナビエ-ストークス方程式によって記述される。ナビエ-ストークス方程式は解析する流体の性質によっていくつかの種類があるが、非圧縮性の流れの場合は式\eqref{NavierStokes}で表される。

\begin{align} \label{NavierStokes} \frac{\mathrm{D}\boldsymbol{v}}{\mathrm{D}t} = -\frac{1}{\rho}\mathrm{grad}p + \frac{\mu}{\rho}\Delta\boldsymbol{v} + \frac{1}{\rho}\mathrm{grad}(\boldsymbol{v}\cdot\mathrm{grad}\mu) + \frac{1}{\rho}\mathrm{rot}(\boldsymbol{v}\times\mathrm{grad}\mu) - \frac{1}{\rho}\boldsymbol{v}\Delta\mu + \boldsymbol{g} \end{align}

ナビエ-ストークス方程式の一般解は求められていない[1]為、流体の解析では主に数値解析を用いて問題を解く。特に流体を数値解析する場合数値流体力学と呼ばれる。

数値流体力学（CFD）

数値流体力学(CFD : Computational Fluid Dynamics)とは流体に関する方程式をコンピュータで計算する事により数値的に流れを求める手法の事である。CFDでは流体の方程式をコンピュータで計算する為、空間や時間と言った連続的な値を離散化して解を求める。離散化された空間ではその一つ一つの離散点に同じ方程式が適用される。

CFDは主に陰的な解法によって計算される。陰的な解法とは計算する方程式の従属変数が他の方程式によって定義されるような場合に反復計算などを行って解を求める解法である。

GPUとGPGPU

GPUとはGraphics Processing Unitのことであり、コンピュータの画面表示を司る装置である。昔はCPUのアクセスできるメモリの一部を用いて画面のイメージを表示していたが、現在のコンピュータでは図形などの描画用のAPIが用意され、PGUが図形などの表示を担当する。そのため、図形表示などの特定の処理を高速に実行できるように画像処理に最適なハードウェアが構築されている。具体的にはPGUは1つの命令で同時に多数のデータを処理するSIMD(Single Instruction Mdltiple Data)と呼ばれる方法で計算を行う。もともとPGUは2Dの描画機能をCPUから切り離しCPUの負荷を低減させる為の物であったが、次第に3Dグラフィックスの描画高速化としてベクトルや行列の演算を行えるよう拡張されていった。ベクトルや行列等の計算が行えるようになると、その計算能力をグラフィックス目的だけでなく汎用的に使用するニーズが生まれてきた。こうしてGPPGUが生まれた。

GPGPU(General-Purpose computing on Graphics Processing Units)とは、PGUを汎用計算に応用する技術の事である。PGUの得意とするベクトルの並列処理をあらゆるデータに対し行えるようにする物で、単純な命令を多量のデータに並列に適応させる事で高速化を図る。なお、GPPGUはあらゆる処理を高速化できるわけではない。例えば一般にGPPGUでは条件分岐は計算速度を損なわせる。これはPGUの演算器に条件分岐を最適化する投機実行や分岐予測と言った機能が無い[2]為である。また、PGUはメモリアクセスの際一定量をまとめて読み込む仕組みの為、アクセスするメモリ位置が連続していない場合は何度もメモリ転送が発生し速度が低下する。したがって、特定の目的をGPPGUで高速化できるかどうかは重要な研究課題である。

図00

SIMD 一回の命令で複数のデータを処理する

CFDとGPPGU

CFDとGPPGUは非常に親和性が高い。これは空間の離散化と言う点にある。CFDでは解析する空間領域を離散化し、各離散点ごとに流体の方程式を計算するが、CPUの場合は2017年現在での並列度は多くて数十なので、個々のコアで全離散点を逐次的に計算するため時間がかかる。一方GPPGUの場合離散点の計算を並列に行うのにSIMD計算が使え、またコア数も数千から数万となるためこの計算時間を大幅に減らすことが出来る。

本研究の概要

本研究では流体の物理シミュレータの高速化手法として、GPPGUを利用した流体シミュレーションプログラムを与えられた格子数やスレッド数などで最適化を施し動的に生成するプログラムを作成した。高速化の手法としてプログラムの動的生成が有用であるか評価する。

2章では本研究で使用したGPPGU環境であるCUDAの基礎的な用語や技術を解説する。

3章では本研究の流体シミュレーション方法として使用した格子ボルツマン法について述べる。

4章では本研究で作成した流体シミュレータのシミュレーション内容について述べる。

5章では本研究の内容についてCPUコードによるLBMの実装からCUDAでの実装、そして本研究で行った高速化手法の解説する。

6章では本研究の結果を示し、評価を行う。

7章では本研究のまとめを述べる。

8章では計測値の生データを示し、本研究で作成したプログラムとその説明を書く。

CUDAについて

CUDA(Compute Unified Device Architecture)とは、NVIDIA社が提供するハードウェア仕様、プログラミングモデル、SDKやAPI等で構成されるGPPGUプログラミングの開発環境の事で、C言語を使用して並列プログラムを作成することが出来る。

CUDAのアーキテクチャ

CUDAにはSMX（Streaming Mdltiprocessor）と呼ばれる装置が多数搭載されている。このSMXにはキャッシュやレジスタ、そして多数のCUDAコアと呼ばれる演算装置等が搭載されている。CUDAコアは計算を行う最小単位と言え、このコアを最大限使用する事がPGUの性能を最大限発揮させる事に繋がる。なおSMX及びSMX内のCUDAコアの数は製品によって異なる。

図00

SMXの構成(Keplerコンピュート・アーキテクチャ・ホワイトペーパーより)

カーネル

CUDAにおいて並列プログラムを作成する場合、並列実行するプログラムをカーネルと呼び、関数として記述される。カーネル内の命令はCUDAコアによって実行される。CUDAでは多数のコアそれぞれが同じカーネルを実行することによって並列処理を行っている。

スレッドとブロック、グリッド、Warp

CUDAではカーネルを管理する為にスレッド階層と呼ばれる階層的な管理方法を提供している。スレッド階層の概要図を図aaaに示す。

スレッドとはカーネルの実行単位である。CUDAでは非常に多数のスレッドを作る事が出来るが、ブロック当たりのスレッド数は上限が存在する。

ブロックはスレッドを束ねたものである。1つのブロックは1つのSMXに割り当てられる。なお1ブロック内に含める事の出来るスレッド数は製品によって異なる。ブロックは最大3次元でスレッドを管理する事が出来る。またブロック内のスレッドは共有メモリを通してデータのやり取りを行う事が出来る。

グリッドはブロックを束ねたものである。グリッドも最大3次元でブロックを管理する事が出来る。ただし、1つの次元で格納する事の出来るブロック数は、製品によって異なるが、上限が存在する。

なおCUDAではグリッドの上位にデバイスと呼ばれる概念が存在する。1つのデバイスは1つのPGUを意味するため、マルチPGU環境では複数のデバイスを持つ事になる。

CUDAではPGU上で動作するそれぞれのスレッドに対しスレッドID（1ブロック内のスレッドの通し番号）とブロックID（1グリッド内のブロックの通し番号）を付与する為、この値を用いてデータの選択を行う事が出来る。これらはカーネル内で組み込み変数として参照できる。

スレッドはWarpと呼ばれる単位ごとに制御される。1Warpは32スレッドであり、Warp内のスレッドはすべて同じ命令を実行する。つまり、Warp内のスレッドはすべて同じカーネルを実行する。また半Warpは16スレッドの事を言う。

図00

スレッド階層

ホストとデバイス

デバイスとは前述のとおり1つのPGUを意味する。一方ホストとはPGUと接続されているCPU（あるいはCPU上で動くカーネルを管理するプログラム）の事を言う。CUDAではホスト側からデバイスへカーネルの実行命令を出すと非同期にデバイスがカーネルの処理を始め、ブロッキングされる事無くホスト側へ処理が戻る。よってホスト側ではデバイス側の処理終了を待つこと無く、他のタスクを実行することが出来る。

メモリ管理

CUDAにはいくつものメモリが存在する。PGUのチップから見たメモリ構造を図aに、スレッド階層から見たメモリ構造を図aaに示す。

図00

チップから見たメモリ構造

図00

スレッド階層から見たメモリ構造

図aaaaを見るとOn ChipとOff Chipとある。これはメモリがPGUのチップのどこに存在しているかを示す。On Chipの場合はチップ内部に、Off Chipの場合はチップ外部に存在している。On Chipメモリの場合は小容量ではあるもののCUDAコアからは高速にアクセスする事が出来る。ただしCPUからの読み書きはできない。一方Off Chipは基本的に大容量であるがCUDAコアからのアクセスはOn Chipメモリよりも低速になる。しかしCPUからのアクセスを行う事が出来る。

レジスタ

レジスタはCUDAコアの最も近くに存在するメモリで非常に高速に読み書きすることが出来る。アクセスする事が出来るのはそのレジスタを割り当てられたスレッドのみとなる。またレジスタを使用するのはカーネル内のローカル変数である。

CUDAではブロック（SMX）内のレジスタを、スレッド（CUDAコア）に分配する。このため、ブロック内のスレッド数が少なければ1つのスレッドが使用する事の出来るレジスタ数は増える。逆にブロック内のスレッド数が多ければ少なくなる。レジスタの数が足りなくなった場合はローカルメモリへデータは退避される。ローカルメモリはレジスタよりも非常に低速なため、レジスタ不足は実行速度の低下を招くので回避する必要がある。その方法はカーネルで使用するローカル変数を減らすか、ブロック内のスレッド数を減らす必要がある。

共有メモリ

共有メモリはOn Chipのメモリで同一ブロック内のスレッドならばどれでも読み書きする事の出来るメモリである。On Chipメモリなので高速にアクセスすることが出来、スレッド間のデータ交換等に使用される。

ローカルメモリ

ローカルメモリはOff Chipのメモリでレジスタや共有メモリよりもアクセス速度は低速になる。ローカルメモリにアクセス出来るのは同じアドレスにつき1つのスレッドまでとなる。ローカルメモリは足りなくなったレジスタを補う為のメモリでCPUやカーネルから明示的にアクセスする事は出来ない。ローカルメモリは図arch-memoryの様に2か所に存在していて、1つはL2キャッシュ、もう一つはグローバルメモリである。最もメモリアクセスが遅いグローバルメモリを使用する前に一段キャッシュ機構を組み込むことでアクセスの高速化を図っている。L2キャッシュが一杯だったりキャッシュミスを起こすとグローバルメモリにアクセスする。

コンスタントメモリ

コンスタントメモリはOn Chipのキャッシュを持つ読み込み専用メモリである。デバイスからの値の書き込みはできないが、CPUから書き込む事が出来る。コンスタントメモリでは値の読み込みはまずOn Chipのキャッシュを通して行われるため、半Warpのすべてのスレッドがコンスタントメモリ内の同じ場所を参照する場合はグローバルメモリを使用する場合に比べて高速に値をロードできる。コンスタントメモリの容量はデバイスによって異なるが、本研究で用いたNVIDIA GeForce TITANでは64KBとなっている。

グローバルメモリ

グローバルメモリはデバイスの中で最も大容量のメモリである。デバイス、ホストの両方で読み書きを行う事が出来る。一方メモリアクセスの速度は最も遅く、実行時間に大きな影響を与える。グローバルメモリアクセスの高速化方法にコアレスアクセスと呼ばれるものがある。これは連続したスレッドがグローバルメモリ内の連続した領域にアクセスする事を言う。CUDAではグローバルメモリアクセスの際連続した領域を一度に読み込むため、連続したスレッドがグローバルメモリの連続した領域にアクセスする場合一度のアクセスで読み込みを行う事が出来る。一方不連続な領域にアクセスする場合、何度もメモリ読み込みを行う為速度が低下する。

グローバルメモリのアクセスを減らすための手段としては、カーネルの最初で値をキャッシュしておく方法が挙げられるがこの場合カーネル内で使用するレジスタの数が増えるため場合によっては一度に並列処理できるスレッド数が減り、結果的に速度が低下する事もある。

図00

コアレスアクセス

格子ボルツマン法

概念

格子ボルツマン法(Lattice Boltzmann Method、以下LBM)はCFDの一種である。これは空間内の気体粒子の分布を離散化された空間で時間発展的に解くことにより、気体の動きを解析する手法の事である。LBMでは空間を有限個の格子に分割し、その格子点の仮想的な粒子の分布を周囲の格子点を用いて計算する。また時間と速度の離散化も行い、離散化された時間をステップと呼ぶ。粒子は1ステップで格子線に沿って別の格子点へちょうど到達する速度で移動する。

LBMでは空間の次元数と静止を含め、向きを考慮した速度の数をそれぞれ$n$と$m$として$n$D$m$Vと言うように表す。例えば、2次元で9つの速度を持つ場合2D9Vとなり下図の様に表すことが出来る。

図00

2D9V

ここで数字は速度を離散化した時の指標である。

格子ボルツマン方程式

LBMは衝突と並進を表す二つの分布関数を計算する事で粒子分布を求める。衝突後の粒子分布を$\overline{f_a}(t,\boldsymbol{r}) $、衝突前の粒子分布を$f_a(t,\boldsymbol{r})$とすると、衝突を表す式は(\ref{collision_step})となる。

\begin{align} \label{collision_step} \overline{f_a}(t,\boldsymbol{r}) = f_a(t,\boldsymbol{r}) + \Omega_a[f_a(t,\boldsymbol{r})] \end{align}

ここで$a$は粒子の持つ速度の方向、$\boldsymbol{r}$ は位置ベクトル、$\Omega_a[f_a(t,\boldsymbol{r})]$ は衝突による粒子分布の変化を表し$\Omega_a$は衝突演算子と呼ばれる。

続いて並進は$\tau$ を時間刻み、$\boldsymbol{c}_a$を$a$方向の並進速度ベクトルとすると式(\ref{streaming_step})となる。

\begin{align} \label{streaming_step} f_a(t+\tau,\boldsymbol{r}+\boldsymbol{c}_a\tau) = \overline{f_a}(t,\boldsymbol{r}) \end{align}

式(\ref{collision_step})と式(\ref{streaming_step})を組み合わせる事によって、衝突と並進による状態変化を考慮した式(\ref{lattice_boltzmann})が求まる。この式を格子Boltzmann方程式と言う。

\begin{align} \label{lattice_boltzmann} f_a(t+\tau,\boldsymbol{r}+\boldsymbol{c}_a\tau) = f_a(t,\boldsymbol{r}) + \Omega_a[f_a(t,\boldsymbol{r})] \end{align}

式(\ref{lattice_boltzmann})の右辺第二項は衝突項と呼ばれる。衝突項の計算方法に本研究では衝突による粒子分布の変化量が一定の割合で減少して平衡状態に近づいていくとするBGKモデルを用いて衝突項を式(\ref{BGK_model})の様に置き換えた。

\begin{align} \label{BGK_model} \Omega_a[f_a(t,\boldsymbol{r})] = -\frac{1}{\phi}[f_a(t,\boldsymbol{r})-f^{(0)}_a(t,\boldsymbol{r})] \end{align}

式(\ref{lattice_boltzmann})と式(\ref{BGK_model})を組み合わせると式(\ref{lattice_BGK_model})が求まる。

\begin{align} \label{lattice_BGK_model} f_a(t+\tau,\boldsymbol{r}+\boldsymbol{c}_a\tau) = f_a(t,\boldsymbol{r})-\frac{1}{\phi}[f_a(t,\boldsymbol{r})-f^{(0)}_a(t,\boldsymbol{r})] \end{align}

式(\ref{lattice_BGK_model})は格子ボルツマン方程式における格子BGKモデルと呼ばれる。ここで$f^{(0)}_a(t,\boldsymbol{r})$は局所平衡分布関数と呼ばれるもので、次のような式である。

\begin{align} \label{peq} f^{(0)}_a(t,\boldsymbol{r}) = w_a\rho[1-\frac{3}{2}(\boldsymbol{u}^2)+3(\boldsymbol{e}_a\cdot\boldsymbol{u})+\frac{9}{2}(\boldsymbol{e}_a\cdot\boldsymbol{u})^2] \end{align}

ここで$\boldsymbol{e}_aと$$w_a$は2D9Vの場合、

表00

$\boldsymbol{e}_aと$$w_a$の値

a	$\boldsymbol{e}_a$	$w_a$
0	(0,0)	4/9
1,2,3,4	(1,0),(0,1),(-1,0),(0,-1)	1/9
5,6,7,8	(1,1),(-1,1),(-1,-1),(1,-1)	1/36

となる。

一方、各格子点の密度$\rho$と流速$\boldsymbol{u}$は次式から求められる。

\begin{equation} \label{density} \rho = \sum_{a}^{}f_a \end{equation}

\begin{equation} \label{velocity} \boldsymbol{u} = \sum_{a}^{}\frac{c\boldsymbol{e}_af_a}{\rho} \end{equation}

ここで$c$は移流速度として$c=\delta x/\delta t$で定義される。$\delta x$、$\delta t$はそれぞれ空間刻み、時間刻みを表す定数である。

式(\ref{lattice_BGK_model})をある格子点の分布関数を求めるよう変形すると

\begin{equation} \label{lattice_BGK_model_2} f_a(t+\tau,\boldsymbol{r}) = f_a(t,\boldsymbol{r}-\boldsymbol{c}_a\tau)-\frac{1}{\phi}[f_a(t,\boldsymbol{r}-\boldsymbol{c}_a\tau)-f^{(0)}_a(t,\boldsymbol{r}-\boldsymbol{c}_a\tau)] \end{equation}

となる。プログラムでは式(\ref{lattice_BGK_model_2})の方を用い、時刻$t$の状態から時刻$t+\tau$の格子点の分布関数を求めている。

シミュレーション内容

本研究では図cavity_flow_simの様な正方キャビティ流れと呼ばれる物理現象をシミュレーションするプログラムを出力する。本章では正方キャビティ流れの説明と共に、シミュレーションに使用した境界条件についても説明する。

図00

キャビティ―流れのシミュレーン内容

正方キャビティ流れ

図cavity_flow{の様な正方形形状の領域を考え、3辺（左右、下部）は固定された固定壁、残る一つ（上部）が一定速度で動く移動壁とする。この時正方形領域内の流体は移動壁に引きずられて動き、右側の壁に衝突して下方向に曲げられる。曲げられた流体はそのまま進み、やがては下部の壁にぶつかり左へ曲げられる。他方、上部左側では流体が右へ移動する為、その不足分を補うために下部から流体が流れ込む。このようにして領域内では最終的に流体が循環するような流れが発生する。この流れの事をキャビティ流れと呼ぶ。

本研究で作成したプログラムはこのキャビティ内流れを計算するカーネルを出力する。

図00

キャビティ―流れの概念図

境界条件

CFDではシミュレーション領域で壁などの流体の運動に影響を与えるものが存在する事がある。その場合は流体の方程式だけでは解を決める事は出来ない為、補助的な条件を設定する事により解を求める。この補助的な条件の事を境界条件と言う。

Bounce-Back境界条件

LBMは式(\ref{lattice_BGK_model_2})より、ある格子点の計算に周囲の格子点の値を利用す。従って、壁面に接する格子点では壁面から流れ込む粒子の分布関数は式(\ref{lattice_BGK_model_2})では求めることはできない。そこでBounce-Back条件を用いて近似する。

Bounce-Back境界条件は粒子を侵入方向から$180^\circ$跳ね返すもので、8番目の方向に進む粒子の場合図bounce_back_graphの様になる。

\begin{equation} \begin{split} f_1(t+\tau,\boldsymbol{r})=f_3(t,\boldsymbol{r}) \\ f_2(t+\tau,\boldsymbol{r})=f_4(t,\boldsymbol{r}) \\ f_5(t+\tau,\boldsymbol{r})=f_7(t,\boldsymbol{r}) \\ f_6(t+\tau,\boldsymbol{r})=f_8(t,\boldsymbol{r}) \\ \end{split} \end{equation}

図00

bounce-back境界条件

本研究では図cavity_flowの左右、下側の三方向の壁面にBounce-Back境界条件を適用している。

ディリクレ境界条件

ディリクレ境界条件とは格子に直接値を与える境界条件である。正方キャビティ流れでは図cavity_clow_2の移動壁を再現する為に領域上部において右向きに一定の流速を与えている。

研究内容

本研究では一つのカーネルに掛かる計算時間を短縮させる為にカーネル内でのメモリアクセスが少なくなるよう調整したカーネルプログラムを出力するプログラムを作成した。

なお、使用した言語はC#で出力されるカーネルプログラムは2次元の正方キャビティ流れをLBMで計算するプログラムである。

本プログラムは次の手順でカーネルプログラムを作成する。

図00

プログラムの流れ

CPUにおけるCFDの実装

CUDAによる高速化を行う前に、基本的なCFDを行うCPUコードの説明を行う。プログラムの流れ図を図cpu_flowに示す。CPUコードはCUDAで並列に計算していた格子毎の計算を全ての格子に対しループを用いて計算している。また、CUDAでは境界条件による値の選択をif文無しに行っていた。これはGPGPUがif文などの条件分岐の処理に対し最適化されておらず、一般に条件分岐は計算速度を損なわせるからである。これに対しCPU側ではそのような制限は無いので格子の種類による計算処理の方向を条件分岐によって決定している。

図00

CPUによるCFDの流れ

CUDAによる高速化

CUDAは非常に多数のスレッドを用いて並列計算を行う事が出来る。CFDでは一般に空間を離散化し、その離散点一つ一つで方程式の計算を行う事により全体の結果を導き出す。よってCUDAでのCFDの計算では、スレッド一つを離散点一つに対応させ全離散点を並列計算する事で高速化を図る。

しかし、CUDAは使用するGPUによって生成できるスレッド数やブロック数が制限される。その為生成できる総スレッド数が総計算点よりも少ない場合が発生する事を考慮し、一つのスレッドで複数の計算点の計算を行うようカーネルを作成する必要がある。

CUDAによるCFDのプログラム

図00

CUDAによるCFDの基本的なプログラム

図base_programはCFDの基本的なプログラムを図示したものである。

計算に使用するデータは現在の状態($t$ステップ)を保持する物(以降、データA)と計算後($t+1$ステップ)の値を格納する物(以降、データB)の二つを用意する。

CUDA側では大量に生成されたスレッドがデータAを用いて$t+1$ステップの状態を計算し、その結果をデータBに保存する。その際、スレッドは複数の格子点をループを用いて計算する。CPU側ではCUDAに計算命令を出した後、計算が終了するまで待機する。計算が終了したらデータAとデータBのポインタを入れ替えて再びCUDAに計算命令を出す。データAとデータBのポインタを入れ替える事でCUDAは$t+1$ステップのデータを用いて状態を計算し始める。その為、結果的に$t+2$ステップの状態を計算できる。この様にデータのポインタを入れ替える事でステップを進めている。

データのやり取りはグローバルメモリを介して行われる。また気体定数等の定数はコンスタントメモリにCUDA初期化時に転送する。

データ構造

出力されるカーネルで使用するデータは一次元配列で管理し、引数として渡される。データは更新前と更新後の二つを用意する。カーネルの引数は格子点の密度、速度、分布関数値を格納するdouble型のポインタdensityとn_densityと、流体セルや壁、流入や流出境界を示すint型のポインタmaskとn_maskを持つ。先頭にn_ を持つ引数は次ステップ用の配列である。

密度、速度、分布関数値を格納する方の一次元配列は、例えば2D9Vの環境で格子が図latticeのように$3\times3$の格子で左上から右下へ0～8までの格子番号が振られているとすると、データの並びは図memoryように先頭に格子番号が0～8番までの密度の値が、その後ろに同じく0～8番までのX方向の風速が並ぶ、このように格子番号順にY方向の風速、Z方向の風速、速度方向0～8番の分布関数値が後ろに続く。

図00

格子構造（格子内の値は格子番号）

図00

データ構造

格子の種類の扱い

通常シミュレートする空間には図lattice_bgk_modelの様に流体の流入口や流出口、壁などが存在する。これらには式(\ref{lattice_BGK_model})を適用できないのでカーネル内で格子の種類により値を変化させる必要がある。そこで格子の種類を整数値の識別番号を用いて表現する事でカーネル内で識別できるようにする。カーネルに渡される格子の種類を格納した配列は図array_mappingの様に空間の左上から右下へマッピングされた1次元配列として渡される。

図00

格子の種類

図00

格子のマッピング

計算領域外の扱い

LBMでは注目している格子の値を周りの格子の値を利用して計算する為、領域の端においては計算を行う事が出来ない。そこで計算領域との境界に境界条件を付加する事で計算を可能にする。

本研究では条件分岐を行う事による計算速度の低下を防ぐため、図boundaryの様に計算領域外との境界をあらかじめ格子として用意し、計算を行うのはその内側（図boundaryにおける青い格子)のみにする事で参照する格子点が領域外かそうでないかの判定を行わずに済むようにしている。

図00

境界格子の構成

計算部分

カーネルプログラムの流れを図program_flowに示す。

カーネルプログラム内ではループを用いて1スレッド多格子点の計算を実現している。ループ内は1つの格子点の計算を行う。まずループ回数を元に計算する格子点のデータ位置を算出する。次に全方向に対する式(\ref{lattice_BGK_model}）、式(\ref{peq})の計算を逐次行う。この時計算格子点の周りの格子点の値を用いて分布関数値を計算するが、その周りの格子点が壁や流入境界、流出境界だった場合を想定して値の選択を行う。値の選択はListingno_branchの様に条件分岐を使わずに行われる。

リスト00

条件分岐を使わない値の選択（例）


結果= (mask == NORMAL) * 分布関数値;
結果+= (mask == BOUNDARY) * 壁(bounc-back)境界の値;
結果+= (mask == INFLOW) * 流入境界の値;
結果+= (mask == OUTFLOW) * 流出境界の値;

ここでNORMAL、BOUNDARY、INFLOW、OUTFLOWはC言語のマスクで、それぞれ流体セル、壁、流入境界、流出境界を表す識別番号に置き換えられる。maskは現在計算に使用している格子点の種類を示す整数型の識別番号で前述の4つのマスクの値のいずれかを持つ。

Listingno_branch_listingはC言語の関係演算子(ここでは同値を意味する==演算子)が真では1を、偽では0を返す事を利用して値の選択を行っている。つまり関係演算子の結果と各条件での値を乗算し、全体を加算する事で関係演算子が真の値のみが結果に保存される。

全方向の分布関数が計算し終えると式(\ref{density})、式(\ref{velocity})を計算し密度と速度を算出、次ステップのデータに保存する。

図00

計算の流れ（2D9Vの場合）

プログラム生成による高速化

カーネルの分割

CUDAではループ回数が異なるスレッドがWarp内に存在する場合は全スレッドで最大のループ回数に合わせてループが実行される。よってWarp内の全てのスレッドが同じ回数ループするのならば問題はないが、異なる回数を実行する場合は図warp_kinituka_の様に無駄なリソースを使うことになる。この場合、Warp内のスレッドが行うループ数を均一化する事でリソースの空きをなくすことが出来る。CFDでは生成可能なスレッド数を分割した格子点数が上回る場合を想定して、1スレッドで多格子点を計算できるようにしている。このためにはカーネル内でループを行うが、格子点数によってはスレッドによってループ回数が異なってくる場合がある。本研究ではこれを避けるために全てのスレッドで同じ回数ループするカーネル（\roop ）と、ループをしないカーネル(\noroop )、余ったスレッドを実行するカーネル(\rem )の最大三つに分割する。分割に際しては、計算格子点数($l$)、上限ブロック数($b_{max}$)、ブロック当たりの上限スレッド数($t_{max}$)を元に実行ブロック数、ブロック当たりの実行スレッド数、ループ回数を計算する。

図00

Warp内の均一化

カーネルのパラメータ計算

ループ回数($r_n$)は式(\ref{roop_calc})で決定される。

下付き文字のnは整数で0の時\roop 、1の時\noroop 、2の時\rem のパラメータを表す。

\begin{equation} \label{roop_calc} r_n = \begin{cases} \lfloor l / (b_{max}* t_{max}) \rfloor & (n = 0)\\ 1 & (n = 1)\\ 1 & (n = 2) \end{cases} \end{equation}

なお$b_{max}$と$t_{max}$はGPUの制限を超えない範囲の任意の整数であり、GPUの能力を最大限生かすためには$t_{max}$は32の倍数である事が望ましい。

実行ブロック数($b_n$)とブロック当たりの実行スレッド数($t_n$)は式(\ref{block_thread_calc})、式(\ref{thread_calc})で決定される。

\begin{equation} \label{block_thread_calc} b_n = \begin{cases} b_{max} & (n = 0) \\ l \boldsymbol (b_{max} * t_{max}) / t_n & (n = 1) \\ 1 & (n = 2) \end{cases} \end{equation}

\begin{equation} \label{thread_calc} t_n = \begin{cases} t_{max} & (n = 0,1) \\ l \boldsymbol (b_{max} * t_{max}) \boldsymbol t_{max} & (n = 2) \end{cases} \end{equation}

カーネル毎の格子データ配列のオフセット($offset_n$)は式(\ref{offset_calc})で決定される。

\begin{equation} \label{offset_calc} offset_n = \begin{cases} 0 & (n = 0) \\ offset_{n-1} + r_{n-1} * t_{n-1} * b_{n-1} & (n = 1,2) \end{cases} \end{equation}

カーネル内の定数展開

各カーネルのパラメータが計算し終わると、パラメータを元にカーネル内の定数を展開する。主に展開するのは配列のインデックス値や方向ベクトルの値である。

Listingnoroop_source、deploy_sourceに本研究で作成した生成プログラムによって出力された定数の展開前と展開後のカーネルコードの一部を示す。

リスト00

展開前（一部抜粋）


index_a = index_out - ex[0] - c_data.x_max_out * ey[0] - c_data.x_max_out * c_data.y_max_out * ez[0];
_vx += ex[0] * c_data.c * pa * boundary;
_vy += ey[0] * c_data.c * pa * boundary;
_vz += ez[0] * c_data.c * pa * boundary;

リスト00

展開後（一部抜粋）


index_a = index_out - (0) - (256) * (0) - (65536) * (0);
_vx += 0 * c_data.c * pa * boundary;
_vy += 0 * c_data.c * pa * boundary;
_vz += 0 * c_data.c * pa * boundary;

ここで展開前のc_dataはコンスタントメモリに格納された格子点の情報を持つ構造体である。またex、ey、ezは方向ベクトルの各要素の配列で、c_data同様コンスタントメモリに格納されている。

このように格子点数等の情報をカーネルに埋め込むことで、コンパイラの最適化を促しさらにメモリアクセスを減らすことで高速化を図る。

生成されるカーネル

格子点数が256×256で使用するブロックとスレッド数を128×128、速度数を9にした時の出力されるカーネルを付録のListingkernel_no_deployment、kernel_deploymentに示す。出力されるプログラムの説明を次に示す。

main.cu: main.cuはlbm_kernel_deployment.cuh及びlbm_kernel_nodeployment.cuhの計算用カーネルを動かすためのメイン関数を持つプログラム。include文の参照を変える事で実行するカーネルを切り替える。
lbm_kernel_deployment.cuhとlbm_kernel_nodeployment.cuh: lbm_kernel_deployment.cuhは定数の展開を行ったカーネルプログラム。; lbm_kernel_nodeployment.cuhは定数展開を行わなかったカーネルプログラム。
lbm_config.cuh、lbm_data.cuh: lbm_config.cuhはカーネル内のマクロを定義し、lbm_data.cuhはプログラム内で使用する構造体を定義している。この二つのファイルは定数の展開の有無で共通である。

評価

測定環境及び測定方法

測定環境を表environmentに、測定条件は表measurement_conditionsに示す。

表00

測定環境

OS	Windows10 64bit
CPU	Intel Core i7-4820K CPU 3.70GHz
GPU	NVIDIA GeForce GTX TITAN
CUDA	Ver 7.5
RAM	32GB

表00

測定条件

次元数	2次元
速度方向数	9
計測ソフト	NVIDIA Visual Profiler Veer7.5

計測は$256\times256$、$512\times512$、$1024\times1024$格子に分割した空間をそれぞれ1000、10000、100000ステップ計算した。またそれぞれの計測は3回行い平均値を結果とする。表中の展開の列は、有りの項目は\deploykernel を、無しの項目は\nodeploykernel の実行結果を意味する。また$256\times256$格子、1000、10000ステップの計算をOpenMPを使用したCPUコードで行った。

結果

それぞれの計算時間の結果を図result_grafuに示す。またCUDAコードにおける$256\times256$格子、1000ステップ計算時のスレッド当たりのカーネル数を図register_threadにOcuppancyの値を図occupancy、グローバルメモリの読み込み回数と書き込み回数をそれぞれ図GLT、図GSTに示す。

図result_grafuから、\deploykernel の方が\nodeploykernel よりも約21%高速に動作する事が確認できた。

また図register_thread、から、定数の展開を行った場合は行わなかった場合に比べ全体的にレジスタの数が減少している。このレジスタ数の減少はカーネルコードにコンスタントメモリに記録されていた値を展開したことによって、コンスタントメモリの値を格納するためのレジスタが削減できたためであると考える。つまりListingnodeploy_sourceにおいてex、ey、ez配列やc_data構造体から値をいったんレジスタに退避していた所を、Listingdeploy_sourceの様に値を直接展開したことによってレジスタに退避する必要がなくなったという事である。

図occupancyが示すOccupancyの向上は、カーネルに必要とされるレジスタ数が減少したことによりSMXが一度に実行できるWarp数が増加した事が原因であると考える。

図GLT、図GSTから展開有りの方が無い場合に比べ読み込み回数が減少している一方、書き込み回数には変化がない。コンスタントメモリはカーネルからは読み込みは出来ても書き込みが出来ない事を考えると、コンスタントメモリの値を展開したことによってコンスタントメモリアクセスが減った事がこれらの数値に影響を与えていると考えられる。またCUDAではOff Chipのメモリアクセスは演算に比べて非常に低速である事から、メモリアクセス回数が減った事が実行時間の削減に寄与していると考えられる。

以上より、カーネルの実行速度が向上したのにはコンスタントメモリを利用していた定数値をカーネル内に展開する事でメモリアクセスの回数が削減されたこと、レジスタ数が減少したことによりSMXで一度に実行できるWarp数が増加したことが原因であると考えられる。

図00

実行時間-ステップ数

図00

スレッド当たりのレジスタ数

図00

Occupancy

図00

グローバルメモリ読み込み回数

図00

グローバルメモリ書き込み回数

まとめ

図result_grafuからCPUコードに比べCUDAコードの実行速度はCUDAコードの高速化を施していない場合であっても約100% の高速化が図れていることが分かる。この事から、CFDの計算をCUDAに行わせるだけでも非常に高速に動作する事が分かった。

本研究ではLBMを用いた流体シミュレータをプログラムで出力する事により、定数をプログラム内に埋め込む事で計算の高速化を図った。結果として約21%の高速化を図ることが出来、事前計算による定数の埋め込みはCUDAの更なる高速化に対して有効であると考える。

本研究では初めにLBMとCUDAを利用した流体シミュレーションカーネルを開発し、その後開発したカーネルコードを元に定数を埋め込むジェネレートプログラムを作成した。このジェネレートプログラムは21%高速化したカーネルを出力できるが、カーネルの内容をプログラム内に埋め込んでいる為非常に汎用性に乏しいプログラムとなってしまった。

今後の課題として、よりカーネルの自由度を高めるためにCUDAのカーネルコードを読み込みその内容を解析して動的に定数を埋め込んだり、あるいは式の事前計算や0との乗算と言った不要な計算を省けるような機構を持つシステムの開発を進めていきたい。

参考文献

David B. Kirk・Wen-mei W. Hwu ,株式会社Bスプラウト (2010) 『CUDAプログラミング実践講座』,株式会社ボーンデジタル
CUDA TOOLKIT DOCUMENTATION』,[online] http://docs.nvidia.com/cuda/index.html
蔦原道久・高田尚樹・片岡武 (1999) 『格子気体法・格子ボルツマン法』,株式会社コロナ社
荒木健・越村俊一・大家隆行 (2005) 『格子ボルツマン法による物体周りの流れの解析』 [online]www.tsunami.civil.tohoku.ac.jp/hokusai3/J/shibu/19/araki.pdf
石川裕・里深信行 (2000) 『格子ボルツマン法による物体周りの流れの数値計算』(第14回数値流体力学シンポジウム資料)[online]www2.nagare.or.jp/jscfd/cfds14/pdf/c08-2.pdf
Massimo Bernaschi・Massimiliano Fatica・Melchionna・Sauro Succi and Efthimios Kaxiras (2009) 『A flexible high-performance Lattice Boltzmann GPU code for the simulations of fluid flows in complex geometries』[online]onlinelibrary.wiley.com/doi/10.1002/cpe.1466/full
P.R. Rinaldi・E.A. Dari・M.J. Venere・A. Clausse (2012) 『A Lattice-Boltzmann solver for 3D fluid simulation on GPU』[online]http://www.sciencedirect.com/science/article/pii/S1569190X1200038X

付録

測定データ

測定データ

プログラム

コードジェネレータ

リスト00

LBMCodeGenerator.cs
リスト00

GeneratorTest.cs

生成されるプログラム（共通プログラム）

リスト00

lbm_config.cuh
リスト00

lbm_data.cuh
リスト00

main.cu

生成されるプログラム（展開無しカーネル）

リスト00

lbm_kernel_nodeployment.cuh

生成されるプログラム（展開有りカーネル）

リスト00

lbm_kernel_deployment.cuh

ジェネレートプログラムの説明

Constantクラス

LBMの局所平衡分布感化数(式(\ref{peq}))の係数wの値を文字列として保持するクラスである。

Directionクラス

方向ベクトルを保持するクラスである。内部はVector3<double>型の配列を持ち、そのサイズは速度方向数分となる。

Invertクラス

Bounce-Back境界条件用の速度番号の反転テーブルを保持するクラスである。内部ではdouble型の配列を持ち、サイズは速度方向数分となる。値にはインデックスの番号と反対方向の速度番号が入っている。

KernelInfoクラス

このクラスは生成されるLBMの計算カーネル(roop_kernel、noroop_kernel、remaining_kernelカーネル)に関する情報を保持するクラスである。

メンバ変数について説明する。

Uint64 roop: カーネル内のループ回数を保持する。
Uint64 blocks: カーネルを実行するブロック数を保持する。
UInt64 thread_par_block: カーネルを実行するブロック当たりのスレッド数を保持する。
UInt64 total_threads: 3つのカーネルで累積されるスレッド数を保持する。データのオフセット値計算に使用。
bool enabled: カーネルを書き出すかどうかのフラグ。ループ回数が1以上かつブロック当たりのスレッド数が1以上の場合trueとなる。

メソッドについて説明する

void dispDebug(): メンバ変数の情報を標準出力に出力する。

LBMConfigクラス

このクラスは格子点数やLBMの速度の数、使用できる最大ブロック数およびブロック当たりの最大スレッド数等の事前計算に必要な情報を保持するクラスである。

メンバ変数を説明する。

Vector3<Uint64> lattice: 各次元の格子点の数を保持する。
Vector3<Uint64> lattice_in: 各次元の計算する格子点の数を保持する
uint direct: 速度の数。
uint speed: 最大速度。
Tuple<String,int>[] masks: カーネル内で格子の種類を特定する為のC言語のマクロ名と識別番号のリスト。
uint max_block: カーネルに使用できる最大ブロック数。
utin max_thread: カーネルに使用できるブロック当たりの最大スレッド数。
LinkedList<KernelInfo> list: 分割するカーネルの情報リスト。
Coefficient w: 局所分布関数(式(\ref{peq}))の係数$w$の値。
Direction e: 方向ベクトル。
Invert inv: Bounce-Back境界条件用の速度番号の反転テーブル。

メソッドについて説明する。

UInt64 getAllLattice(void): 格子点の総数を返す。
void setLattice(uint,uint,uint): 各次元の格子点の数を設定する。引数は初めからX方向、Y方向、Z方向である。
String getLBMDataCode(String): 現在保持しているパラメータでLBMData構造体を初期化するプログラムコードを生成し返すメソッド。第1引数はLBMData構造体の変数名。この変数名を用いて初期化するプログラムコードは生成される。
ulong getOffset(): データ配列における最初の非境界格子までのオフセット値。
void dispLatticeInfo(): 非境界格子の総数を標準出力に出力する。

HostGeneratorクラス

HostGeneratorはCPU側で動作するコード(main.cuファイル)を出力するクラスである。LBMCodeGeneratorを継承する。

HostGeneratorはプログラムのエントリポイントとなるmain関数、速度のデータをファイルに書き込むwriteVelocityData関数、圧力のデータをファイルに書き込むwritePressData関数を出力する。出力はオーバーライドしたgenerate(LBMConfig)メソッドを呼び出すことで行われる。出力中にエラーが発生するとfalseが、完了するとtrueが帰る。

ConfigGeneratorクラス

ConfigGeneratorクラスはプログラムで使用するマクロを定義するlbm_config.cuhファイルを出力するクラスである。generateメソッドの引数にLBMConfigクラスのインスタンスを与える事で格子点数やデータのオフセット値等の定数をマクロとして書き出す。出力はオーバーライドしたgenerate(LBMConfig)メソッドを呼び出すことで行われる。出力中にエラーが発生するとfalseが、完了するとtrueが帰る。

DataGeneratorクラス

DataGeneratorクラスはプログラムで使用する構造体等を定義するlbm_data.cuhファイルを出力するクラスである。生成される構造体は次の通りである。

LBMPoint
LBMResult
LBMConfig
LBMData
LBMPackage

出力はオーバーライドしたgenerate(LBMConfig)メソッドを呼び出すことで行われる。出力中にエラーが発生するとfalseが、完了するとtrueが帰る。

LBMCodeGeneratorクラス

LBMCodeGeneratorクラスはコードを生成しファイルに出力するクラスである。generateメソッドの引数にLBMConfigクラスのインスタンスを渡す事で、LBMConfigクラスの内容に基づいたカーネルの分割、定数の展開を行いファイルとして出力する。戻り値はbool型で生成中にエラーが発生するとfalseが、生成が完了するとtrueが帰る。

メンバ変数としてKernelInfoクラスのリストを保持する。

メソッドについて説明する。

bool generate(LBMConfig): CUDAのプログラムを生成し、書き出すメソッド。内部ではまず引数で与えられたLBMConfigをcalcThreadRoop(LBMConfig)メソッドに与え、分割するカーネルのパラメータリスト(KernelInfoのリスト)を得てをれを引数のlistフィールドに代入する。次にConfigGenerator、DataGenerator、KernelGeneratorDeployment、KernelGeneratorNoDeployment、HostGeneratorクラスのそれぞれに対しgenerateメソッドを引数であるLBMConfigのインスタンスを与え実行し、コードの出力を行う。generateメソッドはコード生成中にエラーが発生するとfalseが帰る為、戻り値がfalseだった時には生成を中断してfalseを呼び出し元に返す。全ての出力が完了するとtrueを呼び出し元に返す。
LinkedList<KernelInfo> calcThreadRoop(LBMConfig): カーネルのパラメータ計算(式(\ref{roop_calc})、式(\ref{block_thread_calc})、式(\ref{thread_calc})、式(\ref{offset_calc}))を実行し、roop_kernel、noroop_kernel、remaining_kernelカーネルの実行ブロック数や実行スレッド数等のパラメータを算出する。そしてそれぞれのカーネルのパラメータを保持したKernelInfoクラスのインスタンスを作成しリストに追加。roop_kernel、noroop_kernel、remaining_kernelカーネル用の3つのKernelInfoクラスをリストに追加し終えたら、そのリストを呼び出し元に返す。

KernelGeneratorDeploymentクラス

KernelGeneratorDeploymentクラスはlbm_kernel_deployment.cuhファイルを出しゅつりょ力するクラスである。このクラスはgenerate関数にLBMConfigのインスタンスを引数に与える事で、格子点数や方向ベクトルの値などをソースコードに展開してファイルを出力する。コードの生成の際にはコードの一部分の生成を担当するメソッドを使用する。これらのメソッドは自分の担当するコードを戻り値として返すため、コード生成はこれらメソッドの戻り値を合成し一つのコードを作り上げる。

このクラスの持つメンバについて説明する。

String argument: roop_kernel、noroop_kernel、remaining_kernelカーネルの引数。

フィールドについて説明する。

virtual bool generate(LBMConfig)

lbm_kernrl_deployment.cuhで定義されるすべての関数を生成し出力するメソッド。内部ではインクルード文、コンスタントメモリの宣言をファイルに書き込んだ後、計算を実行するroop_kernel、noroop_kernel、remaining_kernelカーネル、初期化用のinitKernelカーネル、結果取得用のlbmResultKernelカーネル、CPU側から呼び出され計算の準備やメモリ確保、解放などを行うexchange、all_exchange、mallcData、preparation、initialize、getResultData、execution、allCudaFree関数をそれぞれを生成するメソッドを用いて文字列として取得しファイルに書き込む。

virtual string generate_Kernel(string,string,string,LBMConfig,KernelInfo)

計算を行うroop_kernel、noroop_kernel、remaining_kernelカーネルを書き出すメソッドである。第1引数はカーネル名、第2引数はカーネルの引数、第3引数は戻り値、第4引数は計算条件を保持したLBMConfigクラスのインスタンス、第5引数はカーネルの実行スレッド数などを保持したKernelInfoクラスのインスタンスが与えられる。カーネルの生成は、第1引数、第2引数、第3引数を用いて関数名と引数リストの文字列を作る。その後、カーネル内で使用する変数の宣言と初期化処理を書き出し、generate_thread_roopメソッドを用いて図program_flowの''計算する格子分ループ''のループ内の処理を書き出す。この時、ループ回数が1以下、つまりループしない場合はgenerate_thread_roopメソッドの第2引数に１を与えて実行し、戻り値に格納された処理の文字列を書き出す。もし1より大きい場合は、generat_thread_roopメソッドの内容をループ処理の文字列で囲む。そしてgenerate_thread_roopメソッドの第2引数にループ変数の文字列を与える。こうする事で最終的に生成されるコードは第2引数で与えたループ変数を参照してデータ位置の計算を行う為、結果的に図program_flowの処理を記述する事が出来る。

virtual string generate_thread_roop(LBMConfig,string,ulong,ulong)

図program_flowのループ内の処理を記述するメソッドである。第1引数は計算条件を保持したLBMConfigクラスのインスタンス、第2引数はループ変数の変数名、第3引数はカーネルが計算する格子データまでのオフセット値、第4引数には生成しようとしているカーネルが担当するスレッドの総数が与えられる。このメソッドは図program_flowのデータ位置の計算処理、AからBまでの処理、速度と密度の計算処理、データの保存処理を生成する。AからBまでの処理はgenerate_dist_roopメソッドを速度方向数文ループし各速度方向の分布関数の計算処理を生成する。生成する際、格子点の総数やオフセット値を文字列として埋め込んで生成する。

virtual string generate_dist_roop(LBMConfig,int)

このメソッドは図program_flowのAからB、つまり格子点の分布関数の計算を行う処理を生成するメソッドである。第1引数は計算条件を保持したLBMConfigクラスのインスタンス、第2引数は生成する速度方向の番号である。このメソッドでは第2引数で与えられた生成する速度方向の番号により、方向ベクトルやデータ配列のオフセット値を文字列として生成する。

virtual string generate_result_kernel(LBMConfig)

結果取得用のlbmResultKernelカーネルを生成する。

virtual string generate_getResultData(LBMConfig)

getResultData関数を生成する。

virtual string generate_malloc_data(LBMConfig)

mallocData関数を生成する。

virtual string generate_initialize(LBMConfig)

initialize関数を生成する。

virtual string generate_exchange(LBMConfig)

exchange関数を生成する。

virtual string generate_all_exchange(LBMConfig)

all_exchange関数を生成する。

virtual string generate_execution(LBMConfig)

execution関数を生成する。その際カーネル呼び出し処理のカーネル名はConstantクラスのkernelNamesの値を使用する。

virutal string generate_init_kernel(LBMConfig)

initKernelカーネルの生成を行う。

virtual string generate_allCudaFree(LBMConfig,String)

allCudaFree関数の生成を行う。その際、第2引数で与えられた内容を末尾に付加する。

virtual string generate_preparation(LBMConfig)

preparation関数の生成を行う。その際、次のメソッドを用いて処理を生成する。

virtual string generate_coefficient(Coefficient): 引数から局所平衡分布関数(式(\ref{peq}))の係数wの計算と、コンスタントメモリに転送する処理を生成する。
virtual string generate_direction(Direction): 引数から方向ベクトルの値をコンスタントメモリに転送する処理を生成する。
virtual string generate_invertTable(Invert): 引数からBounce-Back境界条件用の速度方向の反転テーブルをコンスタントメモリに転送する処理を生成する。
virutal string generate_config(LBMConfig): 引数から計算条件の内容を保持したLBMData構造体を初期化し、コンスタントメモリに転送する処理を生成する。

virtual string getFileName()

ファイル名を取得する。このクラスの場合は''lbm_kernel_deployment''が帰る。

このクラスはlbmkernelnodeployment.cuhファイルを出力するクラスである。このクラスはKernelGeneratorDeploymentクラスを継承していて、一部のメソッドをオーバーライドしている。

オーバーライドしているメソッドの変更部分について説明する。

string generatethreadroop(LBMConfig,string,ulong,ulong): データ位置の計算部分において事前計算をせずにコンスタントメモリ内の値を参照し、カーネル実行時に計算させるようになっている。
string generatedistroop(LBMConfig,int): データ位置のオフセット値を事前計算せず、カーネル実行時に計算させるようにしている。
string getFileName(): ファイル名を取得する。このクラスの場合は''lbm_kernel_nodeployment''が帰る。

生成されるプログラムの説明

main.cu

このファイルにはデータの書き込み用関数とエントリポイントとなるmain関数が含まれる。

main.cuの5行目のinclude文は参照するファイルをlbm_kernel_deployment.cuhにすると最適化されたカーネルを、lbm_kernel_nodeployment.cuhにすると最適化されていないカーネルを使用して計算する。

main関数はコマンドライン引数として次の値を取得する。

終了ステップ数
書き込み開始ステップ数
書き込み間隔
空間刻み
時間刻み
代表速度
代表長さ
代表格子点数
分子量
温度[K]
粘性係数
モル気体定数
平均自由行程

main関数は初めにpreparation関数を実行して定数の初期化を行う。

次にinitialize関数を実行しLBMPacage構造体を初期化して、データのポインタを取得する。

初期化が終わると計算を開始する。ループを用いてコマンドライン引数で指定された終了ステップ数まで処理を続ける。

CUDAに計算の実行を指示するのはexecution関数でこの関数が呼ばれるとCUDAは1ステップの計算を開始する。計算が終わるまでブロッキングされる。

計算が終わるとall_exchange関数を呼び出してLBMPackage内のポインタを入れ替えてステップを進める。

もし、現在のステップ数が書き込み開始ステップ数以上で、書き込み間隔との剰余が0ならばデータのファイル書き込みを行う。

その場合getResultData関数を用いてデータを取得する。データはLBMPackage構造体のh_resultメンバに格納される。

LBMPackage構造体のh_resultメンバをwriteVelocityData関数及びwritePressData関数に与えてファイルの書き込みを行う。

ステップが終了ステップまで到達すると、allCudaFree関数を実行し、CUDAの為に生成したメモリなどの解放を行う。

最後にcudaDeviceReset関数を実行しCUDAを完全終了させる。

なお、前述に出てきたpreparation関数、initialize関数、execution関数、all_exchange関数、getResultData関数、allCudaFree関数はいずれも5行目のinclude文で指定したファイル(lbm_kernel_deployment.cuhかlbm_kernel_nodeployment.cuh)の中で定義された関数である。

lbm_config.cuh

このファイルにはプログラム内で使用するマクロが定義されている。マクロの説明を行う。

DIRECTION: 速度方向の数。2D9Vならば9となる。
NORMAL: 格子の種類を表す識別番号。式(\ref{lattice_BGK_model})の計算を行う格子を表す。
INFLOW: 格子の種類を表す識別番号。流入条件を適用させる格子を表す。
OUTFLOW: 格子の種類を表す識別番号。流出条件を適用させる格子を表す。
BOUNDARY: 格子の種類を表す識別番号。Bounce-Back条件を適用させる格子を表す。
RBLOCK: roop_kernelカーネルの実行ブロック数を表す。
RTHREAD: roop_kernelカーネルの実行スレッド数を表す。
NRBLOCK: noroop_kernelカーネルの実行ブロック数を表す。
NRTHREAD: noroop_kernelカーネルの実行スレッド数を表す。
REBLOCK: remaining_kernelカーネルの実行ブロック数を表す。
RETHREAD: remaining_kernelカーネルの実行スレッド数を表す。
LATTICE_X: 格子点のX方向の数。
LATTICE_Y: 格子点のY方向の数。
LATTICE_Z: 格子点のZ方向の数。
LATTICE_IN_X: 境界格子を含まないX方向の格子点数。
LATTICE_IN_Y: 境界格子を含まないY方向の格子点数。
LATTICE_IN_Z: 境界格子を含まないZ方向の格子点数。
DENSITY_OFFSET: データ配列の密度が格納されている位置までのオフセット値。
VX_OFFSET: データ配列のX方向の速度が格納されている位置までのオフセット値。
VY_OFFSET: データ配列のY方向の速度が格納されている位置までのオフセット値。
VZ_OFFSET: データ配列のZ方向の速度が格納されている位置までのオフセット値。
DIST_OFFSET: データ配列の分布関数値が格納されている位置までのオフセット値。

lbm_data.cuh

このファイルにはプログラム内で使用する構造体が定義されている。

LBMResult構造体

データ出力用に1格子点の流速と密度のみを持つ構造体。

double vx,vy,vz: その格子点の流速。
double density}　: その格子点の密度。

LBMConfig構造体

格子点数や気体定数など定数の情報を保持する構造体。

double inflow_a[DIRECTION]: 流入条件で使用するあらかじめ計算された各速度方向の分布関数値を格納する配列。
double tau: 緩和時間。
double c: 移流速度。
double density: 初期密度。
double vx,vy,vz: 初期速度。
double ivx,ivy,ivz: 流入速度。
int max_thread_par_block: 1ブロック当たりのスレッド最大数。
int direct_num: 速度方向数。
int max_speed: 最大速度。
int x,y,z: 各次元の最大数。
int offset: 最初の計算点までのオフセット。
long size_out: 格子点の数。
long size_in: 格子点の数から境界格子の数を引いたもの。
int x_max_out: X方向の格子数。
int y_max_out: Y方向の格子数。
int z_max_out: Z方向の格子数。
int x_max_in: 境界格子を含まないX方向の格子数。
int y_max_in: 境界格子を含まないY方向の格子数。
int z_max_in: 境界格子を含まないZ方向の格子数。

LBMPackage構造体

LBMPackage構造体はCPU側でステップを進めるために行う、データ配列のポインタを保持する構造体である。

LBMResult *h_result,*d_result: h_resultはCPU側のLBMResult配列のポインタを、d_resultはGPU側のLBMResult配列のポインタを意味する。これはGPU側からCPU側へ結果を転送する時に使用する。
double *d_density,*d_n_density: 密度、流速、分布関数の値を保持する配列のポインタ。d_densityは現在の、d_n_density
double *d_mask,*d_n_mask: 格子の種類を表す識別番号を保持する配列のポインタ。d_maskは現在の、d_n_maskは次ステップ用の配列である。

lbm_kernel_nodeployment.cuh

このファイルには最適化を施していないカーネルコードが定義されている。

グローバル変数について説明する。なお先頭に__constant__が付く変数はCUDAのコンスタントメモリに置かれる変数である。

__constant__ LBM::LBMData c_data: 各次元の格子点数や、総格子点数などの定数を持つ変数。
__constant__ double w[DIRECION]: 局所平衡分布関数(式(\ref{peq}))の係数wの値を保持する配列。
__constant__ double ex[DIRECTION]、ey[DIRECTION]、ez[DIRECTION]: 方向ベクトルの値を保持した配列。
__constant__ double invert_table[DIRECTION]: Bounce-Back境界条件による方向の反転用。値は添字の番号と反対方向の番号が格納されている。

各関数について説明する。なお、先頭に__global__ 修飾子が付く関数はGPUで実行できるカーネルプログラムである。

__global__ void roop_kernel(double*,double*,int*,int*)

ループを行い、1スレッド多格子の計算を行うカーネル。図\ref{program_flow}を実装する。第1引数は現ステップのデータ配列、第2引数は次ステップのデータ配列、第3引数は現ステップの格子の種類を表す識別番号の配列、第4引数は次ステップの格子の種類を表す識別番号の配列である。このカーネルはプログラムの動的生成による定数の展開が行われていない。

__global__ void noroop_kernel(double*,double*,int*,int*)

ループを行わないカーネル。図\ref{program_flow}のループを省いた処理を実装する。第1引数は現ステップのデータ配列、第2引数は次ステップのデータ配列、第3引数は現ステップの格子の種類を表す識別番号の配列、第4引数は次ステップの格子の種類を表す識別番号の配列である。このカーネルはプログラムの動的生成による定数の展開が行われていない。

__global__ void remaining_kernel(double*,double*,int*,int*)

実行スレッド数と実行ブロック数の計算の過程で余った格子を計算するカーネル。ループは行わない。第1引数は現ステップのデータ配列、第2引数は次ステップのデータ配列、第3引数は現ステップの格子の種類を表す識別番号の配列、第4引数は次ステップの格子の種類を表す識別番号の配列である。このカーネルはプログラムの動的生成による定数の展開が行われていない。

__global__ void initKernel(double*,int*)

現ステップのデータをc_dataで与えられた初期値で初期化するカーネル。第1引数は初期化する現ステップのデータ配列、第2引数は初期化する現ステップの格子の種類を表す識別番号の配列。

__global__ void lbmResultKernel(double*,LBM::LBMResult*)

第1引数に与えられたデータの速度と密度を分離して第2引数に与えられたLBM::LBMResultの配列にコピーするカーネル。データの取得時に使用する。

void exchange(void**,void**)

第1引数に与えられたポインタと第2引数に与えられたポインタを入れ替える。all_exchange関数で使用される。

void all_exchange(LBM::LBMPackage&)

第1引数で与えられたLBMPackage構造体内のd_densityメンバとd_n_densityメンバ、d_maskメンバとd_n_maskメンバのポインタを入れ替える。CPU側からステップを進めるために呼び出される。

void mallocData(double**,int**)

CUDAのGPU側のメモリ確保用関数cudaMalloc関数を用いて計算に必要なデータ配列と格子の種類を表す識別番号用の配列を確保し、そのポインタの値（GPU側のポインタの値）を第1引数と第2引数の実態に代入する。initialize関数で現ステップと次ステップ用のデータ確保のために使用される。

void preparation(LBM::LBMConfig)

このメソッドでは次の値をコンスタントメモリに転送する。

局所平衡分布関数(式(\ref{peq}))の係数w。
方向ベクトル。
Bounc-Back境界条件用の反転テーブル。
LBMData（シミュレーション環境を表す定数）の値。

void initialize(LBM::LBMPackage&)

このメソッドでは第1引数とmallocData関数を用いて現ステップと、次ステップ用のデータを確保する。その後CUDAでinitKernelカーネルを実行しデータの初期化をCUDAに行わせる。データ取得用のGPU側のメモリを確保し、第1引数のh_resultメンバにポインタを代入する。以上の処理を行い、計算を開始する準備を整える。

void getResultData(LBM::LBMPackage&)

計算の結果を取得する為の関数。一旦CUDAで動いているすべてのスレッドの計算が終了するのを待ち、lbmResultKernelカーネルを実行してデータ取得用のメモリに現ステップから速度と密度を代入させる。その後GPU側のデータ取得用メモリの内容をCPU側のメモリである第1引数のh_resultメンバにコピーする。

void execution(LBM::LBMPackage&)

1ステップの計算をCUDAに実行させるカーネル。 roop_kernel、noroop_kernel、remaining_kernelの3つのカーネルを実行する。その際、生成時に計算されC言語のマクロとして定義されている実行スレッド数と実行ブロック数を指定する。なおroop_kernelの実行スレッド数はRTHREADマクロ、実行ブロック数はRBLOCKマクロに定義されている。noroop_kernelの実行スレッド数はNRTHREADマクロ、実行ブロック数はNRBLOCKマクロに、remaining_kernelの実行スレッド数はRETHREADマクロ、実行ブロック数はREBLOCKマクロに定義されている。

void allCudaFree(LBM::LBMPackage&)

GPU側のメモリをすべて開放する関数。第1引数内のGPU側のメモリが格納されているすべてのメンバに対してCUDAの関数であるcudaFree関数を使用してメモリ解放を行う。

lbm_kernel_deployment.cuh

このファイルには最適化を施したカーネルコードが定義されている。

lbm_kernel_nodeployment.cuhとの違いは、roop_kernel、noroop_kernel、remaining_kernelカーネル内において動的生成によって方向ベクトルや格子数などの定数がソースコード上に展開されている点である。その他のinitKernel、lbmResultKernel、exchange、all_exchange、mallocData、preparation、initialize、getResultData、execution、allCudaFree関数はすべてlbm_kernel_nodeployment.cuhと同じ処理を行う。

[1] ナビエ-ストークス方程式にはミレニアム懸賞問題として100万ドルの懸賞金がかけられている。

[2] PGUはコアをチップ上に大量に入れる一方、CPUに存在する投機実行や分岐予測と言った機能を省きコア数を稼いでいる。

プログラム自動生成とCUDAによる格子ボルツマン法を利用した流体シミュレータの高速化

目次

はじめに

流体と流体力学

数値流体力学（CFD）

GPUとGPGPU

CFDとGPPGU

本研究の概要

CUDAについて

CUDAのアーキテクチャ

カーネル

スレッドとブロック、グリッド、Warp

ホストとデバイス

メモリ管理

レジスタ

共有メモリ

ローカルメモリ

コンスタントメモリ

グローバルメモリ

格子ボルツマン法

概念

格子ボルツマン方程式

シミュレーション内容

正方キャビティ流れ

境界条件

Bounce-Back境界条件

ディリクレ境界条件

研究内容

CPUにおけるCFDの実装

CUDAによる高速化

CUDAによるCFDのプログラム

データ構造

格子の種類の扱い

計算領域外の扱い

計算部分

プログラム生成による高速化

カーネルの分割

カーネルのパラメータ計算

カーネル内の定数展開

生成されるカーネル

評価

測定環境及び測定方法

結果

まとめ

参考文献

付録

測定データ

プログラム

コードジェネレータ

生成されるプログラム（共通プログラム）

生成されるプログラム（展開無しカーネル）

生成されるプログラム（展開有りカーネル）

ジェネレートプログラムの説明

Constantクラス

Directionクラス

Invertクラス

KernelInfoクラス

LBMConfigクラス

HostGeneratorクラス

ConfigGeneratorクラス

DataGeneratorクラス

LBMCodeGeneratorクラス

KernelGeneratorDeploymentクラス

生成されるプログラムの説明

main.cu

lbm_config.cuh

lbm_data.cuh

lbm_kernel_nodeployment.cuh

lbm_kernel_deployment.cuh