RI
STニュースNo．
６４（２０１８）

ビット演算によるCFD（数値流体力学）と等価な高精度流体解析手法
AFlui
dAnal
ysi
sMet
ho
dbybit
wiseo
per
atio
nsforac
hievi
nghi
ghacc
urac
yo
fCFD
高度情報科学技術研究機構
松岡浩
流体解析の時間発展計算をビット演算で超高速に実行する方法として“格子ガス法”がある。
まず、流体が存在する空間中に格子を張り、多数の仮想粒子を格子点上に配置する。仮想粒子は、格子点上でのみ他の仮想粒子と衝突して進行の向きを変えながら格子点間を移動していく。
この時間発展過程をビット演算で計算しながら、ときどき、仮想粒子の数や速度を局所平均してマクロな物理量の空間分布を算出し、流体挙動のスナップショットを得る。このとき、仮想粒子がもつ質量や運動量が衝突の前後で保存されるような粒子衝突を想定する限り、その挙動は、自然界におけるある条件下の流体挙動とかなり似たものになる。しかしながら、仮想粒子の速度が、自然界の流体分子のように連続的な値をとれないため、格子ガス法が導くマクロな挙動は、連続流体を仮定している一般的な数値流体力学（CFD）が導く流体挙動とは多少異なったものになる。２５年も前の話であるが、テシャラ（ChristopherM.Teixeira）は、ある工夫をしてこの“多少の差異”を解消する方法を考案した。これによって、格子ガス法は、CFDと同等な精度をもった流体解析手法となりうる。本稿では、テシャラの方法を筆者なりの解釈で説明する。
ある。従って、激しく変化する乱流も仮想粒

１．はじめにもし、流体シミュレーションの時間発展計子の衝突規則から自己組織化して生成できる算をビット演算で直接行うことができれば、かもしれない。概念が簡単であるため、いろ様々なメリットがある。プログラムを工夫す

いろな分野への応用アイデアの発想を誘発でれば、ひとつの格子点に関係する演算を１きる。そして、非平衡系統計力学など今後のビット幅で行うことができる。このため、一発展が期待される分野と深く関係し、サイエ般的な倍精度計算を行う８バイトレジスタをンスとしての魅力も奥深いと思われる。

使用した場合、１CPUコアでも６４格子点についての並列計算を実行できる。また、各時刻２．連続速度分布が可能な仮想粒子モデル

ステップの物理量を実数で追跡する一般的なファインマンの物理学講義によれば、も数値流体力学（CFD）に比べると必要なメモし、後世にただひとつの知識だけしか残すこリ容量を小さくすることができる。さらに、とができないのであれば、“粒子仮説”を選ぶ

実数演算に比べて回路を構成するトランジスとある。
“粒子”の挙動によって諸現象を説明

タ数が桁違いに少ない論理演算とシフト演算しようとする試みは、非常に広範な分野をカが動作の主体なので、消費電力を少なくするバーできる。ここでは、仮想的な粒子が、いことができる。特に重要なことは、ビット演ろいろな向きにいろいろな速さをもって自由算自体では誤差が生じないので、どんな流体に飛びまわれる物理空間を想定する。
“仮想挙動に対しても安定的に答えをだせることで粒子”について以下の仮定をする。

－17－

RI
STニュースNo．
６４（２０１８）（１）連続速度分布が可能な仮想粒子モデルのの３つを想定すれば、次式が成立する。

仮定①仮想粒子は不生不滅である。
②仮想粒子は、空間中を単独で、いろいろな向きにいろいろな速さ

従って、仮想粒子の挙動に関して、３つので並進移動保存則が成立することを導ける。

することができる。すなわち、仮想粒子がもつ速度の向きと大きさ

・仮想粒子の質量保存則：がとる値には、無限個の可能性があり、連続的な速度分布をもつ。

・仮想粒子の運動量保存則：③個々の仮想粒子は、相加的な性質をもつ“物理量”を担うことができる。その物理量が存在する空間位置は、仮想粒子の移・仮想粒子のエネルギー保存則：動にともなって一緒に移動するが、並進移動の過程においては、その量自体は不変である。

すなわち、「連続速度分布が可能な仮想④仮想粒子は、ときどき、同じ局所的な時粒子モデル」では、上記の手順で得られた空間に存在する他の仮想粒子と“相互作３つの保存則は、CFD（数値流体力学）に用”を起こすことができる。このとき、おける“質量保存則（連続の式）”、“運動量お互いの仮想粒子が担う“相加的な物理保存則（ナビエ･ストークス方程式）”、“エ量”の量を変化させることができるが、ネルギー保存則”と同じ代数構造をもち、相互作用に関係した仮想粒子が担う当該保存量の値も同じになることがわかる。

物理量の総和は、相互作用の前後で不変である。

３．離散速度分布の場合の仮想粒子モデル
第２章で述べた仮想粒子モデルは、仮想粒（２）保存則の成立子が持ちうる速度局所的な時空間領域について、向き と速さは無数の可能性がある連続分布だと仮定しをもつ仮想粒子がひとつの格子点にた。本章では、速度の種類を有限個に限定し度において、速存在する粒子数を、その平均値をた離散分布の場合のモデルを考えてみる。こで表わせば、仮想粒子の“数”にれが、“格子ガス法”のモデルに他ならない。

関する保存則として、次式が成立する。
（１）
仮想粒子モデルによる“静止流体”の記述
平衡状態にある“静止流体”を仮想粒子
ただし、テシャラによる格子ガス法の場モデルで表現することを考える。直感的で合は、４次元面心超立方体格子を用いるのはあるが、多数の仮想粒子が“並進移動”で、と“相互作用”を等方的に繰り返していてルは、４次元空間における勾配ベクトになる。

マクロな視点では見かけ上変化がない状態

上記保存則から、個々の仮想粒子が担う“相加的な物理量”として、“質量量”、“エネルギーを想像すればよい。このとき、格子ガス法”“、運動計算用の格子点においては、速度”－18－

仮想粒子の存在確率分布をもつは、速度RI
STニュースの向きによらず、速さごとに一定でNo．
６４（２０１８）

同じ効果をもたらす「別の複数の異なるある。

速度ベクトル“”をもつ仮想粒子の存在確率のセット」を用いることにな（２）仮 想 粒 子 モ デ ル に よ る“運 動 流 体”のる。すなわち、仮想粒子の存在確率分布記述をもとに算出したマク

我々が現実世界で目にする実在流体の挙ロな相加物理量と同じ値が得られるよう動は千変万化するものである。しかし、こな仮想粒子の存在確率分布

れを上記（１）の静止流体に対する仮定とのセットを作りだすことが目標になる。

関連づけるため、次の仮定をする。
①流体のマクロな速度は、について（３）マクロな速度がの“連続関数”である。

仮想粒子の存在確率分布②従って、十分小さい局所的な時空間領域マクロな速度がを考えれば、当該領域に存在する流体のマクロな速度の流体を記述する
で与えられる流体を記述する仮想粒子の存在確率分布

は、一様かつ定常では、第４章で導出するが、 のあると仮定できる。

関数という形で次式のように書ける。

③この局所領域内に存在する格子ガス法計
算用の格子点における速度をもつ仮想粒子の存在確率分布
状態にある静止流体を速度“は、平衡

この後の計算で、この式からマクロな物理”

量を導きだし、CFDの結果と比較するが、で動く座標系から見た場合の存在確率分その際には、布に等しいと考えられる。の向き（ここで、は、向きの単位ベクトル）を変化させる必要

従って、はない。そこで、を、 の大きさでの右辺に現れる仮想粒子の存在確率分布ある“速さ” のみの関数とみて、をもとに、マクロな相加物理量を算出すあることから、テイラー展開を行い、 にれば、CFDにおける３つの保存則を正確対する依存性を分析する。さらに、後述のに再現できることになる。

計算処理において、「での３次精度まで計しかしながら、仮想粒子がもちうる速算すること」及び「格子の空間構造の対称度が有限個の離散速度モデルの場合は、性に起因する仮想粒子がもちうる速度ベク仮想粒子がもつ速度ベクトルとして、トルの対称性から奇数次数が消えること」“を考慮し、”という値を使うことが一般的には許されない。
“について、 に関”といする０次の項（定数項）と２次の項のみのう速度ベクトルが、使うことが許されて存在を仮定し、現時点で以下の形で記述しいる別の速度ベクトル“”にたまたまておく。

一致すればいいのだが、はいろいろな値をとりうるので、一致は不可能である。従って、数値計算では、「速度ベクトル“”をもつ仮想粒子の存以上の結果を、表式１にまとめる。

在確率」を用いることはできず、それと

－19－

RI
STニュースNo．
６４（２０１８）（４）テシャラによる格子ガス法１６９速度モデト）”、“エネルギー密度”、“運動量カレントルによる考察（応力テンソル）”、“エネルギーカレント”

第２章（２）に示した３つの保存則の式に正しく対応し、CFDにおける３つの保存中に現れる５つのマクロな物理量：則を正確に再現できるはずである。
しかし、「離散速度分布の場合の仮想粒子モデル（すなわち、？＝有限個）」に基づ
く仮想粒子存在確率分布を用いて５つのマクロな物理量を求めた場合は、何かの差異が出てくるであろう。このあと、表式１を５つのマクロな物理量に代入して、CFDとの差異を導出してみる。
ただし、以下の計算では、具体性を増すについて、「連続速度分布が可能な仮想ため「テシャラの格子ガス法１６９速度モデル粒子モデル（すなわち、和の上限値？＝（すなわち、？＝１６９）」の場合を想定する。

∞）」に基づく仮想粒子存在確率分布を用いこのモデルの基本的な仮定を「解説１：テれば、それらは、それぞれ、CFDにおけるシャラの格子ガス法１６９速度モデルにおけ“質 量 密 度”、“運 動 量 密 度（質 量 カ レ ンる仮定」に示す。

－20－

RI
STニュース（５）テシャラの格子ガス法１６９速度モデルにNo．
６４（２０１８）な精度をもちうる数値シミュレーション手よるマクロな物理量の計算法になりうる」ことを確信するためには、テシャラの格子ガス法１６９速度モデルに具体的な式変形手順まで理解しておく必要基づいて、“仮想粒子存在確率分布”がある。本稿では、あえて、詳細を示した。

から、５つのマクロな物理量を計算してみなお、式の表示をできるだけ簡単にするよう。結果を表式２から表式６に示す。

ため、本稿では、以下の簡略化表示を用いこれらの式変形は複雑であるが、「格子ている。

ガス法が、ある工夫をすれば、CFDと等価

－21－

RI
STニュースNo．
６４（２０１８）

ただし、：定数

また、上記以外で、以下の式変形に登場する“D”は、面心超立方体格子の次元の数 でD＝ ４、れる定数である。

－22－

は、表 式 １ に 現

RI
STニュース－23－

No．
６４（２０１８）

RI
STニュースNo．
６４（２０１８）（６）パラメータ間の相互関係

ト演算で超高速に実行することができ、各これまでの変数や定数のパラメータの相時刻ステップのシミュレーション計算にお互関係を図１にまとめて示す。

いては、各格子点に到着した仮想粒子の最終的には、テシャラが導入した調節パラメータ（ビットパターンが衝突規則の入力ビットパ定することで、）をある値に設ターンに一致するかどうかを判定する際のに、エネルギーがの仮想粒子の存在３つの条件を同時に満足できる。このとを隠す“粒子マスク”と非存在（ホール）き、格子ガス法のシミュレーションからを隠す“ホールマスク”を有効（０）また時々刻々得られる各格子点における各速度は無効（１）にする乱数発生の確率の比を、の仮想粒子数から、これを局所時空間で平上記の均することにより、すべてのマクロな物理量を高精度（値として設定すればよい。

以上が、格子ガス法によりビット演算をの３次精度）に算出するこ用いてCFDと同等精度の流体解析を行うとができる。

テシャラの方法である。

また、仮想粒子数の時間発展計算はビッ－24－

RI
STニュース

図１４．仮想粒子存在確率分布

パラメータ間の相互関係の導出

ける仮想粒子の存在状態は、有限個のビッ

最後に、第３章（３）の冒頭で先送りしたト列（“状態ベクトル”）で完全に表現する証明を示す。

ことができる。具体的には、速度ベクトル十分局所的な時空間領域
そこには、単位ベクトル
速さ

No．
６４（２０１８）を考えれば、をもつ仮想粒子が存するときは、の向きにマクロなをもつ一様定常速度応するビットを「１」にし、に対をもつ仮想の流体が粒子が存在しないときは、 に対応する存在すると考えてよいと仮定する。この領域ビットを「０」にすればよい。この有限個内に配置された各格子点において、速度をビットの状態ベクトルは、仮想粒子の存否もつ仮想粒子が存在する確率：

を情報を記述するもので、すべての物理量求める。

（質量、運動量、エネルギー等）は、この状態ベクトルから算出できる。また、各格子（１）格子点における仮想粒子の存在状態の点 が と り う る 状 態 ベ ク ト ル の 種 類 は、２有限個通りの可能性がある。

記述
格子ガス法では、各格子点において、仮想粒子がもちうる速度ベクトルの種類は有（２）仮想粒子の衝突による状態ベクトルの限個である。そして、通常のモデルでは、

遷移「各格子点上では、ひとつの速度ベクトル（エネルギーがで速度の向きがここでは、局所時空間領域）

する一様定常速度に存在の流体中に配置をもてる仮想粒子は高々１個である。（たされた格子点の状態ベクトルをだし、静止粒子の個数は、複数でもよい。
）」

で表す。各格子点では、到着した粒子どうという仮定をする。従って、各格子点におしが相互作用（衝突）を起こし、お互いが－25－

RI
STニュースNo．
６４（２０１８）もつ速度ベクトルを変化させて、近傍の格つりあい”）」と仮定する。すなわち、次式

子点に向けて出発する粒子となる。すなわが成り立つ。

ち、到着粒子の存否情報を記述する状態ベクトルが、出発粒子の存否情報を記述する状態ベクトルへと遷移することになる。この衝突による状態遷移の規則を“衝突規則”

ここで考えている格子ガス法モデルでと呼べば、“衝突規則”は、入力となる状態ベクトル

は、そ の 衝 突 規 則 で、を、出力となる状態ベクトを等しく設定するので、に遷移させるとき「１」、させないル

とき「０」になる２値数：。
さらに、到着粒子の状態ベクトルが、衝で表現できる。ただし、すべての状態遷移の可能性を考えると、この２値数は、（２有限個）有限個×（２） 通り必要である。状態遷移を確率的に行う場合は、と突規則の入力状態ベクトルに一致する

確率は、個々の速度ベクトルについて、到着粒子分布が衝突規則の入力状態に一致の値を確するか否かの確率の積で、次式のように与率的に「１」または「０」に変化させればえられると仮定する。
（“分子混沌の仮定”）よい。
この場合、その平均値をで表現する。［解説１］
テシャラの格子ガス法１６９速度モデルにおける仮定：の５．に述べたとおり、通常の格子ガス法では、衝突規則を、が成立するように設定する。

ベクトル以上のことから、各格子点においては、到着粒子の状態ベクトルの中の、速度ベクトルをもつ仮想粒子の存否情報を表現する
ビットの数値を表わし、「１」または「０」では、それが衝ある。また、「１１＝１、１０≡１、０１＝０、００

突規則のいずれかの入力状態ベクトルに一致したとき、出力状態ベクトルは、それぞれ、状態ここで、へ≡１」であることに注意すれば、の遷移規則が適用され、それが出発粒子の状態ベクトルになる。

がは、とに一致「１」で、それ以外のときは「０」になるこするとき「１」、一致しないとき「０」にな、その平均値をる２値数をとがわかる。なお、

で表現すれば、到着粒子の状態ベクトルが、衝突によっての全ビットが一致したときのみは、の平均値を表す。に遷（４）仮想粒子とホールの存在を隠すマスク処

移 し、そ れ が 出 発 粒 子 に な る 確 率 は、で与えられる。

理の導入
衝突規則を適用する前に、“エネルギー（３）局所平衡状態の仮定の仮想粒子が存在することを隠すマ十分局所的な時空間領域
クロな流体速度がでは、マスク”の局所平衡状態がと“エネルギーの仮想粒子が存在しないこと（ホールが存在すること）成り立ち、「局所平衡状態では、衝突によっを隠すマスク”て互いに移りかわれる仮称粒子の状態ベク値数であり、通常は｢１｣で何の効果も発揮トルの間の遷移はつりあっている。
（“詳細しない。マスク効果を発揮させるときの－26－を導入する。とは２RI
STニュースNo．
６
４（２０１８）

み、値を「０」にする。シミュレーション

で、「参考」に述べた理由から１と計算では、ある一定の確率で「０」の値をの一次結合で表現できる。

とらせる。このとき、「を、「とが１になる確率」が１になる確率」をで表すと、

上記（３）の最終式は、次式になる。

ここで、１次結合の比例定数はのみに依存すると考えられる。
（５）衝突不変量の導出

上式を、以上の仮定から、局所平衡状態では、次式が成立していると見なせる。

について解くと、マクロな速度がの局所時空間内の格子点における、速度をもつ仮想粒子が存在する確率として、次式が得られる。

この両辺をで割って、５．おわりにテシャラ（ChristopherM.Teixeira）は、彼のMITにおける博士論文：“Continuum Limit

ここで、

は、
“
「が１になる確率と「 が１になる確率 」の比”で、」

ofLatt
iceGasFluidDynamics”,MIT,1993。の最後を次のように締めくくっている。

両辺の自然対数をとると、積は和に変換され、次式が得られる。

A necessary prerequisitefor accuratesimulationofanyo
ftheseapplicationsistheremovalo
ft
hedi
sc
r
eet
ne
s
sart
i
f
act
sf
r
omt
hebas
i
cc
o
ns
er
vat
i
o
ne
quat
io
nso
fl
att
i
ce
gas met
ho
ds．I
ti
so
nl
yt
hent
hat t
he、をも

pro
mis
eo
f～１０００t
i
me
sco
mput
at
io
nale
f
f
i
-

つ入力粒子または出発粒子が存在するとき

co
me
s
c
ie
nc
yi
mpr
ove
me
nt over CFD beのみ「１」の値をとり、存在しない場合は

meani
ngf
ul． Wit
ht
hemo
del
sdeve
lope
din「０」であるから、上式は、存在する各粒子

thi
swo
rk，thisste
phasbeenaccomplished．

についてのは、そ れ ぞ れ、速度の和が衝突の前後内容を意訳的に補足解説すれば：

で保存されること、
すなわち、１格子点１ビット幅で時間発展計算を実行が“衝突不変量”になっていることを示しできる格子ガス法は、CFDに比べて、１/１０のている。

メモリサイズで１００倍近い高速計算を実行できる。しかし、その計算モデルは、実際の流（６）速度をもつ仮想粒子が存在する確率：の導出

体粒子がもちうる連続的な速度ベクトルを離散化しているため、計算結果に“物理的に意味のない量（discreetnes
sartifacts）”が付加が“衝突不変量”になるのされてしまう。これを消し去るのがテシャラ

－27－

RI
STニュースNo．
６４（２０１８）の方法であり、これによって、はじめて、なる。

CFDに比べてほぼ１０００倍の効率の良さをも筆者としても、本手法の今後の実用化をめ

つ格子ガス法が重要な役割を果たせるようにざして、努力を継続中である。

－28－