元々檀さんに誘われて2003年にやったシミュレーションなんですが、 卒論で並列関係の勉強をするのに合わせて書き直してみました。 逐次版と並列版がありますが、並列化のライブラリには Phoenixを使っています。 一年前はずいぶん苦労したのに、今回はさくっと書けました。 少しは成長したかな、と思えて嬉しいような。
DLAは、例えば雪の結晶が成長していくようすを説明する理論です。 雪が出来るとき、なぜ氷の固まりができないのか不思議に思ったことはありませんか？ 空気中の水蒸気は、0℃以下になるといきなり固まるのではありません。 たいていの水蒸気は、他の氷とかチリなどの粒にぶつかって初めて氷になります。 一度氷になると、その分子はそこから動くことはできません。 このシミュレーションでは、中心に一個の氷の核を置き、周囲から一個ずつ水蒸気の 分子を近づけていきます。各粒子はランダムな方向に動かし、 氷の核にぶつかるかマスの外に出たら次の粒子に移ります。 氷になった場合はその場で固まって、色つきのマス目になります。 このシミュレーションを、特に拡大して見ると、なぜ雪が氷の固まりにならないのかが 分かります。水蒸気分子は始めにぶつかった所で固まってしまうので、奥のほうの 空洞は埋められずに残ってしまいます。結果、あの雪の「ふわふわ感」が生まれるわけです。
このシミュレーションでポイントになるのは、「粒子の運動」と「粒子を出現させる場所」です。 一応シミュレーションの正確さと実行速度の兼ね合い…となるのですが、 遅い方法が必ずしも正確とは限らないのが難しいところです、 まず粒子の運動ですが、実際には様々な粒子に衝突して「不規則な」運動をします。 この不規則さを計算機上で実現するため、 今回は時間1ステップごとにX軸、Y軸に沿って1マス運動させることにしました。 また、ランダムウォークをしていると領域の外に出てしまうことが起こりますが、 こうした場合、その粒子については追跡を中止するのか、あるいは右端に到達した分子は 反対側の壁から出現するのかという選択肢があります。 また、シミュレートする範囲についても、結晶が小さいうちは結晶からあまりに離れた部分での 粒子の運動はあまり意味が無いため、シミュレートする範囲を限定してしまう方法もあります。 これを表したのが下図です。
この方法はシミュレートが劇的に速くなりますが、後に述べる粒子の発生方法 との兼ね合いで、使えない場合もあります。 次に出現させる場所ですが、元々の仮定では「粒子は一様に低密度で分布している」となっています。 結晶化が進むにつれて周囲を自由運動している粒子は減ってきますが、それは拡散の性質により 周囲から補給されます。 この様子をシミュレートする場合、新しい粒子を
[A] 全平面内に一様に発生させるか
[B] 領域の最外周部(遠くの周上)で発生させるか
[C] 結晶の近くで一様に発生させるか
[D] 結晶の近くの周で発生させるか
という4種類に加えて、外周を四角にするか円にするかで、色々な可能性が考えられます。
[A]「全体一様分布」は、粒子数が領域に比べて少ない時はうまい仮定になっていると思われますが、 粒子数が多くなってくると、既に結晶が出来ている部分に突如粒子が出現することが何度も起こり、 仮定としては正しくないと思われます。また、この仮定は粒子が運動する領域が大きいことが前提となり、 ランダムウォークの時間が長くなり、シミュレーションに時間がかかります。
[B]「最外周部」は、元々の粒子の密度が薄く、外部からの流入が大部分である場合をうまくシミュレートしていると思われますが、 計算に時間がかかります。
[C]「結晶の近くで一様発生」は、ランダムウォークのさせ方として 粒子が動く範囲を、結晶の周囲だけに限定する場合、 シミュレートの始点もその周囲内に限定するのが自然です。 この方法は計算が高速で進みます。
[D]「結晶の近くの周上」は、ランダムウォークの範囲を限定する場合の 自然な方法として考えました。 これも計算は高速です。
前回チャレンジした時には、「発生は遠くの円周上から」「ランダムウォークは全範囲」 ただし、ランダムウォークについて「上/左には行かない」ものと「下/右にはいかないもの」 に分けてあります。これは元にしたソースがそうなっていてそのままにしていたのですが、 結果を見ると何度シミュレートしても結果が類似していて、まずかったようです。 今回は条件としては「通り抜け」(右端に出たら左から出る)を行う関係で 周/領域は全て正方形としましたが、その他については4種類の条件を全て試してみました。
うーん、やっぱなんか変だ。線みたいに伸びすぎ。
実際に結晶が成長する様子が見れるアプレットです。 modeで上に述べたモードを、sourceで「実際に計算する(realtime calculation)」のか、「データを見る(data A / B)」のか選んで startをクリックすると、描画されます。 もっと大きい絵が見たい!ときは下のsideの値を大きく、もっとゆっくり見たい!ときはintervalを長くして下さい。 mode 0 / 1 / 2 / 3 source realtime calulation / view calculated data misc side: pixel / interval: msec / リアルタイム計算で大きい結晶作るのは、けっこう時間かかります。
逐次版も並列版も dla.ccにまとめてあります。(C++です)
#define MODE 1
という部分を変えると、上のA-Dのアルゴリズムを変更できます。(0->A, 1->B, 2->C, 3->D)
#define PHNX
とすると、Phoenix対応になり、複数マシンでの並列計算ができます。 実行方法は、各プロセッサで
$ ./dla [プロセッサ番号] [プロセッサ数]
というコマンドを実行してください。(対故障的機能は使っていません) Phoenixのインストールは このあたりを参考に。
C++版を元に、 即席Java版もかいてみた。C++とJAVAの違いとか、なかなか勉強になります。あ、並列は無しです。 さらに、もっとJavaっぽくオブジェクト使いまくりで書いてみたのが これら。 構成としては、 まず結晶の点のデータを持つクラスとして Point.javaがあり、 このPointをnext()により順に返す DlaDataSource.javaというインターフェースがあります。 このインターフェースの実装として、 DlaCalculator.java(実際にシミュレーションを行う)と DlaDataReader.java(シミュレーション済みのデータをURLConnectionで取得) を書き、 DlaApplet.javaでは引数により二つを使い分けるようになっています。 いや、これだけの内容にインターフェースとか使ってるのは趣味の世界ですが。