Wernerモデルを使った砂丘の生成

はじめに

この記事は、Wernerモデルを使った風紋の生成 の続きです。

前回の記事では、Wernerモデルを使った風紋の生成方法を紹介しました。今回は、Wernerモデルのアルゴリズムを一部変更することで砂丘のパターンを生成します。

今回の記事で使用するコードは下記のリポジトリで公開しています。

使い方

主要なオプションのみ示します。 --help オプションで詳細を確認できます。

• 広い領域をシミュレーションしたい場合にグリッドサイズを大きくしてください
• 時間発展を長く追いたい場合にイテレーション数を大きくしてください
• 砂丘生成モードを利用するには -d オプションをつけてください
• 環境で利用できる砂量を -q オプションで指定してください

$ python simulate.py -h
usage: simulate.py [-h] [-w V] [-g N] [-i N] [-d] [-q Q] [-s N]
options:
  -g N, --grid_size N   grid size (default: 100)
  -i N, --num_iter N    number of iterations (default: 500)
  -d, --dune            apply settings for dune generation
  -q Q, --quantity Q    sand quantity. valid on dune generation (default: 1.0)

使用例

デフォルトで実行すると、砂量 1.0 で計算します。風はx方向（左から右）に吹いています。

python simulate.py -d -g 140

ランダムな高さの表面から三日月形のバルカン砂丘が生成され、風下に向かって流れていきます。クリックすると動画が開きます。

砂量を増やすとどうなるでしょうか。砂量 3.0 で計算してみます。

python simulate.py -d -q 3.0 -g 140

一列に並んだ尾根が形成され、隣り合う尾根と繋がったり離れたりしながら、風下に流れていきます。環境に供給される砂量が少ないとバルカン砂丘が、多いと横断砂丘が生成される観測事実と整合しています。

動作概要

Wernerモデルを使った風紋の生成 に示したアルゴリズムと大まかなフローは同じです。周期境界条件が設定されたグリッドに対して、差分法を適用します。風紋生成と同様に、次の時間ステップのグリッドを求めるために、2種類のルールを適用します。後述するように、ルールの内容が風紋生成とは一部異なっています。

for _ in range(num_iter):
    apply_jumping_rule()
    apply_rolling_rule()

apply_jumping_rule() 関数は風によって砂が飛ばされる過程を模擬しています。以下の現象を考慮します。

• 勾配が大きいと砂の跳躍距離 jump_length が短い。勾配が小さいと砂の跳躍距離が長い
• 勾配が大きいと飛砂量 sand_quantity が多い（上り坂）。勾配が小さいと飛砂量は少ない（下り坂）
• 砂が無い場所（$z < 0$）から砂は飛ばない

apply_rolling_rule() 関数は、砂が崩れる過程を模擬しており、拡散方程式の拡散項に相当する処理です。高さを均一にする効果があります。風紋のパターン生成と同じルールです。

パラメータチューニング

jump_length と sand_quantity の関数形を、

z = \alpha + \beta \ \tanh \left( \gamma \ \frac{\partial z}{ \partial x} \right)

と仮定しています。$\alpha$, $\beta$, $\gamma$ はパラメータです。求める量が2つあるので、合計6個のパラメータがあります。今回は砂丘のパターンが明瞭に生成されるパラメータを見つけてハードコーディングしています。総当たりでパラメータを見つけましたが、数式を解析することでパラメータを見つけることもできると思います。

    def apply_jumping_rule(self) -> None:
        g = self.grid
        for _ in range(int(0.2 * g.x_size * g.y_size)):
            i, j = random.randint(0, g.x_size - 1), random.randint(0, g.y_size - 1)
            if g[i, j] > 0.0:
                jump_length = int(11.2 - 8.0 * np.tanh(8.0 * g.partial_x(i, j)))
                sand_quantity = 0.5 + 0.4 * np.tanh(2.0 * g.partial_x(i, j))
                g[i, j] -= sand_quantity
                g[i + jump_length, j] += sand_quantity

詳細はリポジトリを参照してください。

おわりに

Wernerモデルを使った砂丘パターンの生成方法と、代表的なパラメータにおける計算結果を示しました。簡単なルールを用いて現実の砂丘によく似たパターンが生成できて面白いと思いました。