第12回JHPCNシンポジウム（2020年7月9日）：2020年度採択萌芽課題研究内容紹介（ポスター）［課題番号]EX20202（萌芽課題研究）
リカレントニューラルネットワークによる高解像度流体解析コードの開発
課題代表者： 技術士事務所AIコンピューティングラボ 松岡 浩

目的と概要：設計アイデアの探索段階で「小型流動制御デバイスから流体機械まるごとスケールまでを手軽に同時解析できる高解像度な流体挙動予測ツール」があれば革新的なものづくりに大いに役立つ。
この実現には、①高解像度を得るために配置する多数の格子点に関する“超並列低コスト計算”②様々なレイノルズ数状態の流体挙動を模擬できる“格子流体粘性の広域発現”の２点が特に重要である。
本研究では、格子ガス法流体解析モデルの時間発展計算を、時間展開されたリカレントニューラルネットワークとしてとらえ、さらにこの計算を１格子点につき１ビット幅の超並列演算で実行する。また、各格子点における仮想粒子の到着状況と衝突散乱後の出発状況を観察し、ある条件が整ったときに個々の仮想粒子に個別の加速操作を行うことで、局所平衡的な自己組織化によって発現するマクロな流体粘性を広い範囲で実現できる可能性を示す。

１．C. M. Teixeiraが考案した高精度な格子ガスモデルから出発する。
本研究は、C. M. Teixeiraが考案した４次元面心超立方格子による「３速さ格子ガス法54速度モデル」を出発点にして、リカレント型ニューラルネットワークの概念を参考にしつつ、高解像度流体解析を可能な限り手軽に実現する方法を探るものである。
Teixeiraの54速度モデルに基づく格子流体解析の結果は、ナビエ・ストークス方程式を解く標準的な数値流体力学(非圧縮性流体)の結果と比較した場合、マッハ数に関する３次の精度まで一致することがTeixeiraの論文において証明されている。
（cf. Christopher M. Teixeira, “Continuum Limit of Lattice Gas Fluid Dynamics”, MIT, 1993）（cf. 松岡浩, 菊池範子, “多速さ格子ガス法実用化展開への手がかり”, SENAC, Vol.49, No.4, pp1-15, 2016）（cf. 松岡浩, “ビット演算によるCFDと等価な高精度流体解析手法”, RIST NEWS, No.64, pp17-28, 2018）

連続体流れ場（連続体の流動）

仮想粒子の格子点挙動（格子点上での衝突散乱と格子点間の並進移動）

格子ガス法流体解析では、流体が存在する空間中に規則正しく格子点を配置し、特定の質量と運動量を担う多数の仮想粒子が、①格子点上での「衝突散乱」、②格子点間の「並進移動」を繰り返しながら移動していく様子を疎視化(平均)して、マクロな流体物理量(密度、運動量、速度等)の挙動を得る。

２．格子ガスの衝突散乱過程は多層パーセプトロンの計算に等しい。
赤丸：仮想粒子、青丸：格子点

各格子点に到着した複数の仮想粒子は、粒子の数と運動量を保存する衝突散乱によって互いにその速度を変え、次の目指すべき格子点に向かって出発する。この過程は、例えば、
・入力層＝到着粒子の存否を表す２値(±１)のビットパターン
・出力層＝出発粒子の存否を表す２値(±１)のビットパターン
・重み＝格子点間の移動の可否を表す２値(±１)で確率的に変化・活性化関数＝微小整数をしきい値とする階段関数で構成する“多層パーセプトロン”が各格子点位置に内蔵されていると考え、その計算を行うことに等しい。

ニューラルネット内蔵格子点のイメージ
123456-11-111-1

到着粒子の存否
入力ビットパターン

入力層 中間層 出力層

多層パーセプトロン
内蔵格子点

1-1
1-1-11123456出発粒子の存否
出力ビットパターン

332４1562４156

到着粒子のビットパターン

出発粒子のビットパターン（-1+1+1-1+1-1）（+1-1-1+1-1+1）

格子点における仮想粒子の衝突変化

左図で、入力層は入力値±１をそのまま中間層に伝える。中間層は入力ビットパターンを識別する。例えば、中間層の一番上に位置する赤色
ニューロンにおいて、赤線入力に+1、紫線入力に－1の重みを設定したとすれば、上から紫(-1)→赤(+1)→紫(-1) →赤(+1) →赤(+1) →紫(-1)の重みパターンになる。このとき、入力層への入力ビットパターンがこの重みパターンと一致するときにのみ積和が最大値(6)になる。これを中間層ニューロンの活性化関数をしきい値＝6の階段関数にすることによって識別する。このニューロンの出力は、出力層のうち活性化すべきニューロン1,3,6にだけ接続しておく。出力層のニューロンは、接続されているシナプスの入力重みをすべて+1に設定し、活性化関数をしきい値＝-1の階段関数にする。このことで、いずれかひとつでも入力が+1になれば活性化（出力＝+1）する。実際の適用に当たっては、いろいろな衝突規則が干渉しないように、このような階層を多段に重ねて適用する。
なお、衝突規則を確率的に適用する場合は、入力層のニューロンを増やし、そこに±１の入力を所定の確率でランダムに与える。

3．格子ガスの時間発展過程はリカレントNNの時間展開計算に等しい。
格子ガスの流れ場全体に関する時間発展過程では、ある格子点が内蔵する多層パーセプトロンの出力層から出力された出発粒子が並進移動して、他の格子点が内蔵する多層パーセプトロンの入力層に到着粒子として入力されることが繰り返される。従って、流れ場全体の時間発展のからくりは、“多層リカレントニューラルネットワーク（RNN）”の動作で表現することができ、時々刻々の過渡変化は、RNNの時間展開計算で求めることができる。
RNNの時間展開（2次元の場合のイメージ）
異なる時刻におけるニューロンどうしが並進移動のルートでシナプス結合されていく。

並進移動と時間推移を可能にする１時刻ステップ
先の異なる格子点間をつなぐｼﾅﾌﾟｽ結合

到着粒子のビットパターン

衝突散乱の計算を行う多層ﾊﾟｰｾﾌﾟﾄﾛﾝ

時刻ステップｔ+1
衝突散乱

時刻ステップｔ
衝突散乱

並進移動

出発粒子のビットパターン

出発粒子
到着粒子
疎視化データ等のスナップショット出力
到着粒子

出発粒子

4．多数の格子点に関する“超並列低コスト計算”を展望する。(個人的見解)
筆者が考える「多数の格子点に関する“超並列低コスト化”」を実現するための道筋を下図に示す。
格子ガス法流体解析の時間発展計算過程は、“多層パーセプトロン”による計算の組合せに他ならない。さらに、各ニューロンの入出力、重み、階段型活性化関数のしきい値をすべて微小整数で設定することが可能である。このため、高解像度計算を行う際の超並列格子点計算を実現するに当たっては、まずは、ニューラルネット計算で高性能を発揮できるAI用高速CPU等の利用が合理的である。専用CPUの選択もありうる。
他方、多層パーセプトロンを構成する個々のニューロン計算を分析すると、前述のとおり２種類のニューロンのみが必要であり、この計算は簡易な論理演算で代行できる。この結果、各格子点ごと各速度ごとに仮想粒子の存否情報を１ビット幅で計算できるので、ベクトル型計算機(例：SX-ACE)なら1CPU(4コア)でも１時刻における演算命令の発信で、理想的には64bits/word✕256words/命令✕4コア＝65536個の格子点に関する演算を一括実行できる。本研究では、この道筋を追求してきた。
なお、次の段階では、確率的に行うビット演算を“ゲート型量子計算機”で瞬時に実行できる可能性も探る。
格子ガス法
流体解析の時間発展
計算過程
4次元FCHC格子に基づく
54速度格子ガスモデル（C. M. Teixeiraが考案したモデルで非圧縮性NS方程式の解法と３次の精度で等価）

多層パーセ
プトロンによる計算の組合せ

微小整数
ニューラル
ネット計算
モデル

AI用高速CPU向き、または専用CPU

１格子点１ビット幅
計算モデル

ベクトル型CPU向き

検討
中ゲート型
量子計算
モデル

ゲート型
量子計算機(常温)の利用可能性

5．全ての機能を含めたリカレントニューラルネットワークを構成する。
格子ガス法流体計算における“衝突散乱”と“並進移動”の計算は、すべてが同じ条件下で実行してよいわけではない。流体中に不浸透な静止した固体壁が存在する場合、通常、“粘着条件”としてバウンズバック条件が課せられる。この場合、到着粒子の衝突散乱計算によって得られた出発粒子(緩和粒子)の存否計算の結果は捨てられ、到着粒子の速度を反転させた存否分布を出発粒子の存否分布とする（下図の緑線）。
計測データ等の外部環境条件入力

並進移動と時間推移を可能にする１時刻ステップ
先の異なる格子点間をつなぐｼﾅﾌﾟｽ結合
到着粒子のビットパターン

衝突散乱の計算を行う多層ﾊﾟｰｾﾌﾟﾄﾛﾝ

緩和粒子のビットパターン

個別粒子への達人操作の計算を行う多層ﾊﾟｰｾﾌﾟﾄﾛﾝ

１時刻ステップ前の出発粒子の存否情報をリカレントする！
疎視化データ等のスナップショット出力

出発粒子のビットパターン

さらに、実世界の流速計測情報をリアルタイムでフィードバックするような場合、外部環境条件に合うように個別粒子を巧みに加速するなど達人的な操作を行い、いわゆる“データ同化” (計測融合)を実現した
い場合もある。また、後述する粘性の広域制御のため、１時刻ステップ前の出発粒子の存否情報を用いる場合もある。
これらをまとめて左図を得る。

6．幅広い流体粘性を発現させる個別粒子への達人操作を導入する。
流体粘性は、結局のところ、ファンデルワールス力等の複雑な分子レベルの機構から発現するため、格子ガスモデルのスケールで意図的な粘性発現機構を考えても“大胆な近似”にしかならないであろう。他方、粘性の効果によって流体挙動が影響を受けるので、流体挙動に応じた粘性発現モデルは意外に近似度が高い実用的なものになるかもしれない。
一般に、格子ガスモデルのスケールでは、個々の格子点における衝突散乱の前後において、到着粒子の頻度分布がどのような出発粒子の頻度分布を生み出すか？が格子流体の粘性を決定すると考えられている。
実際、格子流体が発現する粘性は、格子点間隔を一定にした場合、衝突規則の詳細に依存して変化する。
しかし、通常の衝突規則は、時刻ｔに同じ格子点に到着した各速度をもつ仮想粒子の存否パターンのみを入力とするものであり、他の要因によっては左右されない。
他方、流体粘性は、流体速度の時間相関に関わるものであるため、本研究では、衝突散乱後得られた緩和粒子の存否情報(前頁の図の左から４番目の空色矢印)を意図的に変化させる操作を行うか否かを決定する情報として、時刻ステップｔにおける到着粒子の存否情報(前頁の図の緑矢印)以外に、時刻ステップt－1における出発粒子の存否情報(前頁の図の赤矢印)を用いることにした。すなわち、１時刻ステップ前の出発粒子の存否情報も「“個別粒子への達人操作”を行う多層パーセプトロン」へのリカレント入力とする。
具体的な“個別粒子への達人操作”としては、“衝突散乱過程”のあと「ある向き(順向き)に出発するはずだった仮想粒子(緩和粒子)を、ある条件が整ったときだけ確率的に逆向きに出発させる」という“達人操作”をすべての向きについて同じ確率で行う。ただし、これは、あるひとつの粒子を意図的に逆向きに加速したことを意味するので、局所的な時空間平均で運動量保存則が破れないように、今回の研究では「順向きの加速と逆向きの加速は必ず交互に行う」ものとした。この操作は、次頁のNNで実行できる。

7．格子流体粘性を幅広く発現できる多層パーセプトロンを考える。
「個別粒子への達人操作の計算を行う多層パーセプトロン」に「格子流体粘性の広域発現操作」を行わせる場合の設定例
達人操作の適用確率で+1の値をランダムに与える

中間層のニューロンは、上から順に後述する「順向き局所主流巻込操作」,「逆向き 局所主流巻込操作」,「順向き近傍粒子引合操作」,「逆向き近傍粒子引合操作」を行うべき条件を識別する。

入力層

中間層

出力層

≧7で発火

≧0で発火

順向き出発粒子の存否出力±１

逆向き緩和粒子の存否入力±１

≧7で発火

≧0で発火

逆向き出発粒子の存否出力±１

順向き到着粒子の存否入力±１≧7で発火逆向き到着粒子の存否入力±１≧7で発火順向き緩和粒子の存否入力±１

１時刻ステップ前の順向き出発粒子の存否入力±１
１時刻ステップ前の逆向き出発粒子の存否入力±１

左図で、入力層は入力値±１をそのまま中間層に伝える。中間層は入力ビットパターンを識別する。赤線は+1、青線は－1の重み設定を意味する。例えば、中間層の一番上に位置するニューロンにおいては、その７つの入力線は、上から赤(+1)→青(-1)→赤(+1)→赤(+1) →青(-1) →赤(+1) →青(-1)の重みパターンになる。このとき、入力層への入力ビットパターンが中間層の重みパターンと一致するときにのみ積和が最大値(７)になる。この状態を中間層ニューロンの活性化関数をしきい値＝7の階段関数にすることによって識別する。このニューロンの出力は、出力層のうち活性化すべきニューロンである上のニューロンにだけ接続しておく。出力層のニューロンは、接続されているシナプスの入力重みをすべて+1に設定し、活性化関数をしきい値＝0の階段関数にする。このことで、いずれかひとつでも入力が+1になれば活性化（出力＝+1）する。
なお、これらのニューロン計算も、簡易な論理演算で代行可能である。

8．円柱後流の過渡変化計算により粘性変化の兆候を観察する。
東北大学サイバーサイエンスセンターのベクトル型スーパーコンピュータSX-ACEを利用し、１CPU(４コア)による40分程度の小規模計算で求められる円柱後流の３次元過渡変化計算の結果（約530万格子点で12800時刻ステップ、円柱軸に垂直なある平面上の流体挙動）を次ページ以降に示す。
３次元空間中に多数の格子点をＸＹＺの各方向に並べ、直方体形状の格子点配列を作る。各格子点は、その内部に４次元目の座標としてR＝0,1,2,3の位置を識別できる自由度をもつとする。これが下図であり、４次元面心超立方体格子を３次元空間へ投影した姿である。今回の数値シミュレーションでは、３次元縮退格子として、Ｘ方向に864個、Ｙ方向に384個、Ｚ方向に16個の格子点を配置した。

流体入口

Z

流体出口

流れZ

ここでは、Rのとれる値をR=0,1,2,3 に限定した“４重３次元格子”を使う．

XR

ZYX
Y

４次元空間
X

Y

[過渡変化シミュレーションの条件]
シミュレーション計算を開始する時刻ステップ0の時点で、各格子点には、そこに存在できる仮想粒子の最大数の20％の数の仮想粒子をランダムな向きに配置する。この結果、疎視化して得られるマクロな流速はゼロであり、流体は、直方体形状の中で静止している。 次に、時刻ステップ１の時点から、＋Ｘ向きの速度をもつ仮想粒子をＸ＝0の位置から注入していく。すると、時刻ステップが進むにつれて、流体全体が＋Ｘ向きのマクロな速度をもつようになる。このとき、＋Ｘ側の先にある直方体出口においては、出口直前に存在する格子点上の仮想粒子配置を、出口直後に存在する格子点の仮想粒子配置にコピーして、出口におけるマクロな流速の勾配がゼロになるという境界条件を近似的に実現した。また、±Ｙ方向と±Ｚ方向には、周期的境界条件を適用した。
この流れの中の入り口に近い位置に、“Ｚ方向の中心軸をもつ無限大の長さの円柱”を置き、その後流に生じる流体挙動を計算した。
次頁以降に、あるZ断面上の疎視化された運動量分布の過渡変化動画の中からスナップショット画像を抜粋表示する。

9．流体粘性を変える“局所主流巻込操作(仮称)”の適用効果を観る。
「局所主流巻込操作（仮称）」を確率的に適用してみる。
マクロな速度がゼロで静止している流体が存在する場所にある格子点上では、移動することが可能なすべての向きから等方的に仮想粒子が到着し、衝突散乱後、すべて向きに等方的に仮想粒子が出発している。もし、流体がある向き（順向き）にマクロな速度をもっていれば、順向き速度成分が正の速度をもった仮想粒子の方が、逆向き速度成分が正の速度をもった仮想粒子よりも、頻繁に到着し、かつ、頻繁に出発している。個々の格子点において、一直線上で相対する前後２つの向きへの仮想粒子の移動を考える。ある瞬間に順向きの到着粒子と出発粒子がともに存在し、逆向きの到着粒子と出発粒子がともに存在しない場合、局所的な主流が瞬間的に順向きであると考える。このとき、格子点から逆向きに出発しようとした仮想粒子（緩和粒子）の一部について、主流に巻き込まれて速度を反転させ順向きの格子点に向かって並進していくように操作する。これを、ここでは、“局所主流巻込操作（仮称）”と呼ぶことにする。
巻込操作適用率0%

巻込操作適用率3%

巻込操作適用率7%(LGA06234-50)(LGA06082-50)(LGA06227-50)

巻込操作適用率13%(LGA06225-50)

巻込操作適用率27%(LGA06228-50)“局所主流巻込操作（仮称）”の適用確率を、0%､3%、7%、13％、27%と変化させた場合、円柱後流の過渡変化挙動がかなり変化する。ここでは、挙動が定常化した時点におけるある瞬間（12800時
刻ステップ目）の疎視化された運動量を示す。同一格子点幅で計算時間は同じ。

10．流体粘性を変える“近傍粒子引合操作(仮称)”の適用効果を観る。
「近傍粒子引合操作（仮称）」を確率的に適用してみる。
個々の格子点において、一直線上で相対する前後２つの向きへの仮想粒子の移動を考える。ある瞬間に順向きの出発粒子と逆向きの到着粒子がともに存在し、逆向きの出発粒子と順向きの到着粒子がともに存在しない場合、局所的には、その瞬間、前方の方が後方よりも多くの仮想粒子が存在していたと考える。このとき、格子点から逆向きに出発しようとした仮想粒子（緩和粒子）の一部について、前方に多く存在する粒子からの粒子間引力が勝ることによって速度を反転させ順向きの格子点に向かって並進していくように操作する。これを、ここでは、“近傍粒子引合操作（仮称）”と呼ぶことにする。
引合操作適用率0.1%

引合操作適用率0.5%(LGA060615-50)

引合操作適用率2%(LGA060613-50)(LGA060610-50)

引合操作適用率7%

引合操作適用率20%(LGA06065-50)(LGA06066-50)“近傍粒子引合操作（仮称）”の適用確率を、0.1%､0.5%、2%、7％、20%と変化させた場合、死水域の面積が変化する。また、このとき、円柱が流体から受けている抗力の相対値は、それぞれ、5709, 4131,3207, 2511, 2638であり、死水域の減少とともに抗力が大きく減少した。

11．まとめ

１．これまでの「RNNによる高解像度流体解析コードの開発」について：（１）高解像度を得るために配置する多数の格子点に関する“超並列低コスト計算”
スライド４に示したとおり、“微小整数ニューラルネット計算モデル”と“１格子点１ビット幅計算モデル”ついては、それぞれが、今後の計算機性能の向上によりさらなる“超並列低コスト計算”を実現するのに有効であろう。

（２）様々なレイノルズ数状態の流体挙動を模擬できる“格子流体粘性の広域発現”
現時点では定性的な評価にとどまるが、円柱後流の過渡変化に関する数値実験から、『時刻ｔにおける衝突散乱の結果各方向に出発しようとする“緩和粒子”に対して、当該格子点に時刻ｔに到着した同方向の仮想粒子と時刻ｔ－１に出発した同方向の仮想粒子の有無に関する情報だけを観察し、ある条件が整った際にある確率で速度を反転させて出発させてやる』という“マクスウェルの魔物”にも似た操作（←ここでは“達人操作”と呼ぶ。）を行うと、様々な流体挙動を自己組織化により生成できることがわかった。このとき、円柱の配置体系、入口流入速度、格子点間隔等はすべて同じであるので、この“達人操作”によって格子流体が発現する粘性が変化したものと考えられる。本手法は、様々なレイノルズ数における流体挙動を少ない格子点数による短時間計算でも模擬できるツールに発展させたいと考えている。

２．今後の研究計画について：
上記の“達人操作”に関する定量的理論的考察を進めるとともに、スライド５に示した全体像を完成させるため、“ニューラルネットワークの学習機能”を活用し、以下のようなデータ同化の課題にとりくむ。
①人工的に発現させた格子流体粘性を実在流体の粘性に一致させる計測融合学習②限られた時空間でしか得られない不完全な計測情報から格子流体挙動を実世界に一致させる学習【謝辞】本研究を進めるにあたり、東北大学サイバーサイエンスセンターのスタッフ及びベクトル型スーパーコンピュータSX-ACEの卓越した性能に多大なお世話になった。ここに感謝の意を表する。