可逆計算で流体シミュレーション

[可逆論理で状態可逆分割セル・オートマトンのシミュレーション](https://qiita.com/tadashi9e/items/0e74caa94d8c81835e1f)に続き、別のセル・オートマトンを実装してみました。

## ブロックセル・オートマトン

分割セルオートマトンとしては、マーゴラス近傍を使ったものもあります。

マーゴラス近傍というのは、ステップ毎に境界位置が切り替わるもので、具体的には二次元格子の縦線・横線それぞれを奇数番目の線と偶数番目の線に分けたとして、縦線横線共に奇数番目の線で区切ったときブロックと偶数番目の線で区切ったときのブロックを交互に一つのセルとみなすものです。

ブロックを構成する要素は縦２個×横２個で合計４個です。マーゴラス近傍で面白いのは、奇偶の切り替えによってセル内の４つの要素のセル中心からの向きが反転することです。これを利用すると、単純なルールで粒子の動きを表現することができます。

例えば、一方に向かって動き続ける粒子を表現するルールは、セル中心から見て反対位置に粒子を動かすことで実現できます。

マーゴラス近傍を使ったブロックセル・オートマトンとしてはビリヤードボールモデルが有名です。ビリヤードボールモデルは、可逆計算を表現できるモデルであり、ビリヤードボールモデル上でトフォリゲートなどを含むブール回路を実現することもできます。

## ビリヤードモデル

ビリヤードボールモデルのルールは例えば、以下のようなものです。

> □ □　→　□ □
> □ □　　　□ □

空なら次も空

> □ □　→　□ ■
> ■ □　　　□ □

ボールが来れば、そのまま対角方向に進める

> □ ■　→　■ □
> ■ □　　　□ ■

ボールが正面衝突したら、互いに９０度向きを替えて弾きあう

> ■ ■　→　■ ■
> □ ■　　　□ ■

ボールが３つ衝突したら、来た向きと同じ方向に戻る

> □ ■　→　□ ■
> □ ■　　　□ ■

ボールが２つ９０度の角度で衝突したら、来た向きと同じ方向に戻る

バリエーションは更に列挙できますが、上記を回転したものとして全てのルールを与えることができます。

### 実装

実装方法は状態可逆分割セル・オートマトンとほぼ同じです。
一つのセルは一つの CPU であると見なし、セル内の４つの要素に対応する４つのレジスタを設けます。奇偶の境界の変化は、CPU 間の通信および CPU 内部の場所の入れ替えという考え方で実装しました。通信によってブロック内の要素の情報を受け取ったあとは、４入力制御ノットゲートを用いたパターンマッチによって次のステップの状態を決定します。

## HPP モデル

実際に上記のようなモデルのシミュレーションを行うと若干の違和感があります。

ボール（■）の動きが運動量保存則を満たしていないのです。

ボール（■）が運動量保存則を満たす形でルールを改変すると

> □ □　→　□ □
> □ □　　　□ □

空なら次も空

> □ □　→　□ ■
> ■ □　　　□ □

ボールが来れば、そのまま対角方向に進める

> □ ■　→　■ □
> ■ □　　　□ ■

ボールが正面衝突したら、互いに９０度向きを替えて弾きあう

> ■ ■　→　■ □
> □ ■　　　■ ■

ボールが３つ衝突したら、正面衝突した２つは元の方向に弾きあい、残る一つはそのまま進む。

> ■ □　→　□ ■
> ■ □　　　□ ■

ボールが２つ９０度の角度で衝突したら、弾きあって進む
（ボールに個性はないので、２つのボールがそのまますり抜けて動き続けているようにも見える）

というものになります。これは格子ガスオートマトンの一種（原型）、HPP モデルです。

気体ガスオートマトン法は流体シミュレーションへのセル・オートマトンの応用です。

普通にセル上の粒子の動きを表示するだけだとただのセル・オートマトンにしか見えません。

![hpp_cells64x64.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/412308/4218e566-ed9a-9d1f-ec8f-0d995e459eb2.gif)

物理的な解釈を考慮して、気体の速度（セル上の場所が粒子の運動方向を表すので適当な領域についてそれらの平均をとる）を表示したり、気体の密度（これも適当な領域について粒子数を数える）を表示したりすると、確かに流体っぽく見えます。

![hpp_flow64x64.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/412308/214ec3f7-6f02-e1a6-06c1-2973488a205f.gif)

赤色：高圧領域
水色：低圧領域
白い線：領域毎の流れの向き

画像サイズの関係で長時間動くところをお見せできませんが、さらに時間を進めて見ていると初期位置と（トーラス上の）反対位置との間で粗密が入れ替わって振動が続くのが見えます。

## FHP モデル

HPP モデルだと、偽の保存量があったり、非等方性が見えたりするのですが、三角格子を使うと等方性を満足できて偽の保存量もある程度取り除くことができます。三角格子を使ったモデルは HPP モデルと呼ばれます。

ということは知っていたのですが、実際に FHP モデルのシミュレーションをやったことはありませんでした（もっと正直に言うと HPP モデルのシミュレーションも初めて）。

三角格子はシミュレーションも描画も面倒臭そうで敬遠していたのですが、せっかくなのでやってみました。

FHP モデルもいくつかあって、

- FHP-I モデル：２体正面衝突と、対称な３体の衝突の散乱を考慮する
- FHP-II モデル：正面衝突以外の２体衝突、静止粒子も考慮する
- FHP-III モデル：３体から５体までの全ての衝突を考慮する

の順にパターンが複雑になります。

当然、複雑にするにはそれなりの理由があって、FHP-I モデルだと流体シミュレーションとしては問題があるのですが、今回はお試しなので FHP-I モデルで実装しました。

結論から言うと、あまり芳しい結果とはなりませんでした。

- おそらく衝突則が不十分。森北出版の「セルオートマトン法」に載っていた FHP-I モデルの解説を頼りに衝突則を書き下したのですが、見た感じ６方向の並進運動が目立って見えていて散乱があまり起こっていないように見えます。
- トーラス形状のフィールドがあるだけ、というのは単純すぎる。流れの中に物体を置く、といった問題設定が必要。

というあたりが問題であろうと思います。

![fhp_cells.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/412308/f5831ad3-88a8-057f-f309-328db62fdb6f.gif)

![fhp_flow.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/412308/1c885886-713c-f3ec-5556-b42d9cac44f6.gif)


## まとめ

量子計算の勉強をしてみても量子計算シミュレータで扱えるのはごく少数の量子ビットに過ぎないので「こんな単純な操作しかできないものが有用なのだろうか？」という疑念を抑えることができませんでした。

そこで、「可逆計算で〜」と題した記事では、一旦「量子」を諦めて「古典」もしくは「純粋状態」のみの扱いに的を絞る代わりに扱えるビット数の制約を外してみる、ということをやってみました。

- [可逆論理での箱玉系シミュレーション](https://qiita.com/tadashi9e/items/3f1a3adde974524b0d69)
- [可逆論理で状態可逆分割セル・オートマトンシミュレーション](https://qiita.com/tadashi9e/items/0e74caa94d8c81835e1f)

可逆計算も、実装アイディアを思いついてはゴミ情報消去ができないことに気づいて放棄、ということを繰り返しており、結局 CA の実装に終始しています。「これだけ面倒くさい制約のある計算方法が使い物になるのだろうか？」とやや悲観的な気分でした。

しかし今回、何故か流体シミュレーションに行き着いてしまいました。
流体などという、本来ランダムに動く粒子の動きによって生まれる現象を可逆計算の上で動く可逆 CA の上で、曲がりなりにもシミュレーションできるというのはなんとも不思議です。

もちろん計算理論だのの理屈の上ではできて当たり前ではありますが「量子計算＝役に立たないもの」という思い込みは若干薄らいだ気がします。