格子ボルツマン法の局所平衡分布関数を自動で展開するプログラム

格子ボルツマン法(Lattice Boltzmann Method, LBM)の局所平衡分布関数(local equilibrium distribution)を自動で展開するプログラム．

数値計算のループを実行する際に，展開された（ベクトル，行列を使わない）式を使うことで計算の高速化が期待できるかも？

格子ボルツマン法の局所平衡分布関数は下式で与えられる．下式は[Timm Krüger, et. al., The Lattice Boltzmann Method, Springer, 2017.]による．

この式に，決められた $w_i, c_{i\alpha}, c_s$ を代入して，整理された式を出力するプログラムを示す．

$\displaystyle{ f^{\mathrm{eq}}_{i} = w_{i} \rho \left( 1+ \frac{ c_{i\alpha} u_{\alpha} }{ c_{s}^{2} } + \frac{ u_{\alpha} u_{\beta} (c_{i\alpha} c_{i\beta} - c_{s}^{2} \delta_{\alpha \beta} ) }{ 2 c_{s}^{4} } \right) \tag{3.54} }$

記号の簡単な説明として，

$w_i$ は粒子の速度ベクトルの大きさ要素
$c_{i \alpha}$ は粒子の速度ベクトルの方向要素
$c_s$ は音速
$\delta_{\alpha \beta}$ はクロネッカーのデルタ．

添え字については，

$i$ は粒子の速度ベクトルの番号
$\alpha, \beta$ は速度方向の座標の番号（0,1,2がそれぞれx,y,zに対応）．

pythonプログラムのソースコードを以下に示す．

#!/usr/bin/env python
# coding: utf-8

# In[1]:


import sympy as sym


# In[2]:


rho = sym.symbols("\\rho")
alpha, beta = sym.symbols("\\alpha \\beta")
cs = sym.sqrt(sym.Rational(1, 3))
u = sym.IndexedBase("u")
c = sym.IndexedBase("c")

w_D1Q3 = sym.Array([sym.Rational(2, 3), sym.Rational(1, 6), sym.Rational(1, 6),])
c_D1Q3 = sym.Array([[+0], [+1], [-1]])

w_D2Q9 = sym.Array(
    [sym.Rational(4, 9)] + [sym.Rational(1, 9)] * 4 + [sym.Rational(1, 36)] * 4
)
c_D2Q9 = sym.Array(
    [
        [+0, +0],  # 移動なし
        [+1, +0],  # +x方向
        [+0, +1],  # +y方向
        [-1, +0],  # -x方向
        [+0, -1],  # -y方向
        [+1, +1],  # +x, +y方向
        [-1, +1],  # -x, +y方向
        [-1, -1],  # -x, -y方向
        [+1, -1],  # +x, -y方向
    ]
)


w_D3Q15 = sym.Array(
    [sym.Rational(2, 9)] + [sym.Rational(1, 9)] * 6 + [sym.Rational(1, 72)] * 8
)
c_D3Q15 = sym.Array(
    [
        [+0, +0, +0],  # 移動なし
        [+1, +0, +0],  # +x方向
        [-1, +0, +0],  # -x方向
        [+0, +1, +0],  # +y方向
        [+0, -1, +0],  # -y方向
        [+0, +0, +1],  # +z方向
        [+0, +0, -1],  # -z方向
        [+1, +1, +1],  # +x, +y, +z方向
        [-1, -1, -1],  # 以下同様
        [+1, +1, -1],
        [-1, -1, +1],
        [+1, -1, +1],
        [-1, +1, -1],
        [-1, +1, +1],
        [+1, -1, -1],
    ]
)


# In[3]:


def f_eq(i, w, c):
    Q = c.shape[1]
    _term0 = 1
    _term1 = sym.Sum(c[i, alpha] * u[alpha] / (cs ** 2), (alpha, 0, Q - 1)).doit()
    _term2 = sym.Sum(
        sym.Sum(
            u[alpha]
            * u[beta]
            * (c[i, alpha] * c[i, beta] - (cs ** 2) * sym.KroneckerDelta(alpha, beta))
            / (2 * cs ** 4),
            (alpha, 0, Q - 1),
        ),
        (beta, 0, Q - 1),
    ).doit()
    return sym.simplify(rho * w[i] * (_term0 + _term1 + _term2))
    pass


# In[4]:


f_eq(5, w_D2Q9, c_D2Q9)

jupyter labで実行時の出力結果：

$\displaystyle{ \displaystyle \frac{\rho \left(3 {u}_{0}^{2} + 9 {u}_{0} {u}_{1} + 3 {u}_{0} + 3 {u}_{1}^{2} + 3 {u}_{1} + 1\right)}{36} }$

式の $\rho$ は密度， $u_0, u_1, u_2$ はそれぞれ，x, y, z方向の速度を表す．

この結果は2次元の9速度モデル(D2Q9モデル)の「5番目の式」： $f^\mathrm{eq}_5$ を表し， [Timm Krüger, et. al., The Lattice Boltzmann Method, Springer, 2017.]の式(3.65)に含まれる $f^\mathrm{eq}_5$ と同じ意味を表す．

f_eq(1, w_D1Q3, c_D1Q3)

を実行時の結果：

$\displaystyle{ \displaystyle \frac{\rho \left(3 {u}_{0}^{2} + 3 {u}_{0} + 1\right)}{6} }$

この結果は1次元の3速度モデル(D1Q3モデル)の「1番目の式」: $f^\mathrm{eq}_1$ を表し， [Timm Krüger, et. al., The Lattice Boltzmann Method, Springer, 2017.]の式(3.64)に含まれる $f^\mathrm{eq}_1$ と同じ意味を表す．

雑感等

音楽，数学，語学，その他に関するメモを記す．

格子ボルツマン法の局所平衡分布関数を自動で展開するプログラム