##はじめに
　前回までは、1次元の移流方程式について考えてきました。今回は、多次元への拡張として2次元の移流方程式について数値計算を行っていきます。今回は、Gauss Hill問題と呼ばれる、ガウス分布が一定速度で移動する様子をアニメーションで描画してみようと思います。

##シリーズ構成
####1. スカラー移流方程式
　1.1 輸送速度が正の場合
　1.2 輸送速度の符号が不明の時の線形問題
　1.3 数値流束を利用した方法
　1.4 輸送速度が未知量であるとき(Burgers方程式)
　1.5 多次元への拡張(2次元スカラー移流方程式)
　1.6 実践的な計算法(TVD方程式)
####2. 移流方程式の時間積分法
　2.1 時間陽解法について
　2.2 時間陰解法について(Crank-Nicholson法、近似LU分解)
####3. 2階偏微分方程式
　3.1 1次元熱伝導方程式の数値計算法
　3.2 楕円方程式の数値計算法
####4. 圧縮性流れの数値計算法
　4.1 オイラー方程式の数値計算法
　4.2 MacCormack法による数値計算法
　4.3 TVD法による数値計算法

##支配方程式
　2次元スカラー移流方程式は、

\frac{\partial{q}}{\partial{t}}+c\frac{\partial{q}}{\partial{x}}+d\frac{\partial{q}}{\partial{y}}=0

ここで、$c$および$d$が正の実数であるときを考えます。一次精度風上法を利用すると、

q^{n+1}_{j,k}=q^{n}_{j,k}-c\Delta{t}\frac{q^{n}_{j,k}-q^{n}_{j-1,k}}{\Delta{x}}-d\Delta{t}\frac{q^{n}_{j,k}-q^{n}_{j,k-1}}{\Delta{y}}

と離散化できます、今回は$c=d=1.0$、$\Delta{x}=\Delta{y}=1.0[m]$として計算を行いましょう。

　では、実際にプログラムを実行してみましょう。

from mpl_toolkits.mplot3d import Axes3D 
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.animation as animation

c = 1.0
d = 1.0

dt = 0.05
dx = 0.1
dy = 0.1

jmax = 41
kmax = 41
nmax = 50

ims = []

def init(dx, dy, jmax, kmax):
    x = np.linspace(0, dx * (jmax-1), jmax)
    y = np.linspace(0, dy * (kmax-1), kmax)
    q = init_q(x, y, jmax, kmax)
        
    return (x, y, q)

def init_q(x, y, jmax, kmax):
    def gausfunc(x, y, mu, V):
        detV = np.linalg.det(V)
        invV = np.linalg.inv(V)
        A = 1.0 / (2*np.pi*np.sqrt(detV))
        dX = (np.array([x, y]) - mu[:, None, None]).transpose(1, 2, 0)
        return A * np.exp(-0.5 * dX[:,:, None] @ invV[None, None] @ dX[:,:,:, None])

    mu = np.array([dx * jmax / 4, dy * kmax / 4])
    V = np.array([[dx * jmax / 40, 0], [0, dy * kmax / 40]])
    X, Y = np.meshgrid(x, y)
    
    q = gausfunc(X, Y, mu, V)[:, :, 0, 0]
    
    return q

def do_computing(x, y, q, dt, dx, dy, nmax, interval = 2):
    X, Y = np.meshgrid(x, y)

    fig = plt.figure(figsize=(10, 7), dpi=100) # グラフのサイズ
    plt.rcParams["font.size"] = 10 # グラフの文字サイズ
    
    # 初期分布の可視化
    ax1 = fig.add_subplot(1, 1, 1, projection='3d')
    ax1.set_xlabel('x')
    ax1.set_ylabel('y')

    for n in range(1, nmax + 1):    
        qold = q.copy()
        for j in range(1, jmax - 1):
            for k in range(1, kmax - 1):
                式(2)
                q[j][k] = qold[j][k] - dt / dx * (
                    c * (qold[j][k] - qold[j-1][k]) + d * (qold[j][k] - qold[j][k-1])
                )

        im = ax1.plot_wireframe(X,Y,q, color = "black", rstride = 1, cstride = 1, linewidth = 0.5)
        ims.append([im])

    ani = animation.ArtistAnimation(fig, ims)
    plt.show()

x, y, q = init(dx, dy, jmax, kmax)
do_computing(x, y, q, dt, dx, dy, nmax)

実行すると、

ガウス分布が移動していく様子が見られました。

##参考文献
藤井孝蔵、立川智章、Pythonで学ぶ流体力学の数値計算法、2020/10/20