pythonを用いた2相格子ボルツマン法の実装 

#格子ボルツマン法について
研究で使う機会があり, 苦しかったので情報を共有できたらと思います。
格子ボルツマン法(Lattice boltzmann method 以下LBMと略記)はボルツマン方程式を速度, 時間, 座標に関して離散化して解くCFDの一種です。
以下が簡単なシミュレーション方法になります。
![Picture2.png](https://qiita-image-store.s3.amazonaws.com/0/227484/d334423c-539e-739c-652e-d76023f66ec5.png)

境界条件を入れる場合はstreamingステップの後に入ります。上図のようにアルゴリズムはstreamingとcollisionを繰り返すだけであり, シンプルです。

<!--
#ハマったポイント
###独特な計算手法
格子気体モデルから派生した手法のため, 差分法を使った一般的なCFDと比べると, 無次元化,
境界条件, 偏微分の方法が大きく異なり, 学習に時間がかかりました。

###発展途上の方法
LBMは2000年ごろから発達してきた手法のため, 今でも様々なモデルが考案されています。
そのため,実装にどのモデルを使うべきかの選定が大変でした。 
--> 

# おすすめの参考資料 
[](上記のLBMの特徴から最近出版され,手法の内容が網羅された本を買うことをおすすめします。)

[この本](https://www.amazon.com/Lattice-Boltzmann-Method-Principles-Practice-ebook/dp/B01MRH7CZW)がLBMの基礎方程式の導出から,無次元化, initialization, 境界条件, NS with force, 2相系あたりまで式の証明と引用文献が多数ありとても良くまとまっていると思います。一番お世話になっています。

以下その他資料まとめ。

無次元化
http://www.timm-krueger.de/downloads/Krueger_Edmonton_scaling.pdf

偏微分
https://arxiv.org/pdf/0710.5470.pdf

境界条件(等速or等圧or等密度)
http://phelafel.technion.ac.il/~drorden/project/ZouHe.pdf

境界条件(Straight_Boundaries)
http://lbmworkshop.com/wp-content/uploads/2011/08/Straight_boundaries.pdf

# 2相系のコード
ある流体がより大きい粘性を持つ流体に注入されたときに現れるパターン
(Viscous fingering)のシミュレーション。(どちらも簡単のためニュートン流体でやってます)
非圧縮性流体の方程式と界面力ありのナビエストークス方程式, Cahn–Hilliard 方程式を解いています。
境界条件は上下とブロックはhalfway bounce back, 左が等速壁, 右はopen boundaryです。
参考文献 https://www.sciencedirect.com/science/article/pii/S0377025715002037?via%3Dihub
![ezgif.com-video-to-gif-compressor.gif](https://qiita-image-store.s3.amazonaws.com/0/227484/56b86094-3c74-96c4-2f41-c9a52ed3f263.gif)

以下メインのコードになります。2相系での実装のポイントは注入口と流出口の境界条件において ２つの速度分布関数f, gに対してそれぞれ別々の計算手法を用いることだと思います。
ブロックの作成部分とブロックのbounce backのコードは[こちら](https://github.com/ebinan92/Fingering_dynamics/tree/master/lattice_boltzmann
)にあります。
上のシミュレーション(390×390)でi7-8700でまわして40分ぐらいかかります。

```python
import numpy as np
import sympy
import matplotlib.pyplot as plt
import copy
import cv2
from create_block import Createblock
from bounce_back import Bounce_back
# for server
#plt.switch_backend('agg')
import math
from multiprocessing import Pool
from scipy.ndimage.morphology import binary_fill_holes
import matplotlib.animation as animation

import math
import time

H = 400  # lattice dimensions
W = 420
MAX_T = 26000
psi_wall = 0.0  # wettability on block and wall
Pe = 150  # Peclet number
C_W = 5.0 * (10 ** (-5)) / W  # conversion width
Ca = 7.33 * 10 ** (-3)  # Capillary number
M = 20.0  # Eta non_newtonian / Eta newtonian
R_Nu = 10 ** (-6)  # physical kinematic viscosity of newtonian
tau = 1 / (3.0 - math.sqrt(3))  # relaxation time
rho0 = 1.0  # non-dimensional pressure
n_non = 1.0  # rho0 power-law parameter
Eta_n = 0.023  # Eta newtonian
R_sigma = 0.045  # physical interfacial tension
C_rho = 1.0 * 10 ** 3  # conversion pressure
v1 = (tau - 0.5) / 3  # non-dimensional kinematic viscosity of newtonian
C_t = v1 / R_Nu * (C_W ** 2)  # conversion time step
# x_array = np.arange(1.0, 1.7, 0.01)* 100
sigma = R_sigma * (C_t ** 2) / (C_rho * (C_W ** 3))  # interfacial tension
u0 = Ca * sigma / (rho0 * v1)  # inlet velocity
xi = 2.0  # interface thickness
kappa = 0.75 * sigma * xi  # interfacial tension
a = - 2.0 * kappa / (xi ** 2)
gamma = u0 * W / ((-a * Pe) * (tau - 0.5))
x = sympy.symbols('x')
print("u0:{}".format(u0))


class Compute:
    def __init__(self, mask):
        self.mask = mask # ブロックのない座標(計算する座標)
        self.e = np.array([np.array([0, 0]) for i in range(9)]) # lattice speed
        self.w = np.array([0.0 for i in range(9)]) # weight 
        for i in range(9):
            if i == 0:
                self.w[i] = 4 / 9
                self.e[i][0] = 0
                self.e[i][1] = 0
            elif i == 1:
                self.w[i] = 1 / 9
                self.e[i][0] = 1
                self.e[i][1] = 0
            elif i == 2:
                self.w[i] = 1 / 9
                self.e[i][0] = 0
                self.e[i][1] = 1
            elif i == 3:
                self.w[i] = 1 / 9
                self.e[i][0] = -1
                self.e[i][1] = 0
            elif i == 4:
                self.w[i] = 1 / 9
                self.e[i][0] = 0
                self.e[i][1] = -1
            elif i == 5:
                self.w[i] = 1 / 36
                self.e[i][0] = 1
                self.e[i][1] = 1
            elif i == 6:
                self.w[i] = 1 / 36
                self.e[i][0] = -1
                self.e[i][1] = 1
            elif i == 7:
                self.w[i] = 1 / 36
                self.e[i][0] = -1
                self.e[i][1] = -1
            elif i == 8:
                self.w[i] = 1 / 36
                self.e[i][0] = 1
                self.e[i][1] = -1
            print(self.e[i], self.w[i])
        self.psi = np.full((H, W), -1.0) # 秩序変数
        self.psi[:, :10] = 1.0 # 安定のためあらかじめ左辺に注入しておく
        self.block_mask = np.logical_not(mask)
        self.psi[self.block_mask] = psi_wall # ブロックの秩序変数を設定
        self.left_wall = np.full((H, 1), 1.0) # 左の注入壁
        self.right_wall = np.full((H, 1), -1.0) # 右の流出壁
        self.gamma = gamma
        self.top_bottom_wall = np.full((1, W + 2), psi_wall) # top and bottom wall
        self.nabla_psix = np.zeros((H, W)) # 秩序変数のxの偏微分
        self.nabla_psiy = np.zeros((H, W)) # 秩序変数のyの偏微分
        self.nabla_psi2 = np.zeros((H, W)) # 秩序変数のラプラシアン
        self.rho = np.ones((H, W))[mask] * rho0  # macroscopic density
        self.ux = np.zeros((H, W))[mask] # 速度のx成分
        self.uy = np.zeros((H, W))[mask] # 速度のy成分
        self.p = np.zeros((H, W))[mask] # 圧力
        self.f = np.array([np.zeros((H, W)) for i in range(9)]) # 速度分布関数
        self.g = np.array([np.zeros((H, W)) for i in range(9)]) # 速度分布関数(2)
        self.feq = np.array([np.zeros((H, W))[mask] for i in range(9)]) # 局所平衡解(for f)
        self.geq = np.array([np.zeros((H, W))[mask] for i in range(9)]) # 局所平衡解(for g)
        self.mu = np.zeros((H, W))[mask] # 化学ポテンシャル
        self.mix_tau = np.zeros((H, W))[mask] # 2つ目の緩和時間
        self.F = np.array([np.zeros((H, W))[mask] for i in range(9)]) # force
        self.nabla_psix = self.getNabla_psix()
        self.nabla_psiy = self.getNabla_psiy()
        self.nabla_psi2 = self.getNabla_psi2()
        mu = self.getMu()
        # self.uy = self.getUy()
        self.p = self.getP()
        self.mix_tau = self.getMix_tau()
        for i in range(9):
            self.f[i][self.mask] = self.getfeq(i)
            self.g[i][self.mask] = self.getgeq(i)

    def getP(self):
        return 1 / 3 * self.rho + self.psi[self.mask] * self.mu

    def getUx(self):
        temp = np.zeros((H, W))[self.mask]
        for i in range(9):
            temp += self.f[i][self.mask] * self.e[i][0]
        ux = (temp + self.mu * self.nabla_psix[self.mask] / 2) / self.rho
        return ux
        # print("ux:{}, uy:{}".format(self.ux.mean(), self.uy.mean()))

    def getUy(self):
        temp = np.zeros((H, W))[self.mask]
        for i in range(9):
            temp += self.f[i][self.mask] * self.e[i][1]
        uy = (temp + self.mu * self.nabla_psiy[self.mask] / 2) / self.rho
        return uy

    def getMu(self):
        mu = a * self.psi[self.mask] * (1.0 - np.power(self.psi[self.mask], 2)) - kappa * self.nabla_psi2[self.mask]
        return mu

    def getMu_plain(self):
        mu = a * self.psi * (1.0 - np.power(self.psi, 2)) - kappa * self.nabla_psi2
        return mu

    def getRho(self):
        return np.sum(self.f, axis=0)[self.mask]  # macroscopic density
        # print("rho:{}".format(self.rho.mean()))

    def getA0(self):
        a0 = (self.rho - 3.0 * (1.0 - self.w[0]) * self.p) / self.w[0]
        return a0

    def getA1_8(self):
        a1_8 = 3 * self.p
        return a1_8

    def getB0(self):
        b0 = (self.psi[self.mask] - 3.0 * (1.0 - self.w[0]) * self.gamma * self.mu) / self.w[0]
        return b0

    def getB1_8(self):
        b1_8 = 3 * self.gamma * self.mu
        return b1_8

    def getfeq(self, n):
        if n == 0:
            feq = self.w[n] * (self.getA0() + self.rho * (
                    3 * (self.e[n][0] * self.ux + self.e[n][1] * self.uy) + 4.5 * (
                    self.e[n][0] * self.ux + self.e[n][1] * self.uy) ** 2 - 1.5 * (self.ux ** 2 + self.uy ** 2)))
        else:
            feq = self.w[n] * (self.getA1_8() + self.rho * (
                    3 * (self.e[n][0] * self.ux + self.e[n][1] * self.uy) + 4.5 * (
                    self.e[n][0] * self.ux + self.e[n][1] * self.uy) ** 2 - 1.5 * (self.ux ** 2 + self.uy ** 2)))
        # print("feq{}:{}".format(n, feq.mean()))
        return feq

    def getgeq(self, n):
        if n == 0:
            geq = self.w[n] * (self.getB0() + self.psi[self.mask] * (
                    3 * (self.e[n][0] * self.ux + self.e[n][1] * self.uy) + 4.5 * (
                    self.e[n][0] * self.ux + self.e[n][1] * self.uy) ** 2 - 1.5 * (self.ux ** 2 + self.uy ** 2)))
        else:
            geq = self.w[n] * (self.getB1_8() + self.psi[self.mask] * (
                    3 * (self.e[n][0] * self.ux + self.e[n][1] * self.uy) + 4.5 * (
                    self.e[n][0] * self.ux + self.e[n][1] * self.uy) ** 2 - 1.5 * (self.ux ** 2 + self.uy ** 2)))
        return geq

    def getLarge_F(self, n):
        f = self.mu * self.w[n] * (1 - 1 / (2 * self.mix_tau)) \
            * (((self.e[n][0] - self.ux) * 3 + self.e[n][0] * (self.e[n][0] * self.ux + self.e[n][1] * self.uy) * 9)
               * self.nabla_psix[self.mask] + ((self.e[n][1] - self.uy) * 3 + self.e[n][1] * (
                        self.e[n][0] * self.ux + self.e[n][1] * self.uy) * 9) * self.nabla_psiy[self.mask])
        return f

    def getMix_tau(self):
        v2 = v1 * M
        mix_v = np.divide(2 * v1 * v2, (v1 * (1.0 - self.psi[self.mask]) + v2 * (1.0 + self.psi[self.mask])))
        mix_tau = 3 * mix_v + 0.5
        return mix_tau

    def udpatePsi(self):
        self.psi[self.mask] = np.sum(self.g, axis=0)[self.mask]
        self.psi[self.block_mask] = psi_wall

    def getNabla_psix(self):
        f = np.zeros((H, W))
        psi_with_block = copy.deepcopy(self.psi)
        psi_with_block[self.block_mask] = psi_wall
        temp = np.hstack((self.left_wall, np.hstack((psi_with_block, self.right_wall))))
        temp = np.vstack((self.top_bottom_wall, np.vstack((temp, self.top_bottom_wall))))
        for i in range(9):
            if i == 1:
                f += 4 * np.roll(temp, -1, axis=1)[1:-1, 1:-1]
            elif i == 3:
                f += -4 * np.roll(temp, 1, axis=1)[1:-1, 1:-1]
            elif i == 5:
                f += np.roll(np.roll(temp, -1, axis=1), -1, axis=0)[1:-1, 1:-1]
            elif i == 6:
                f += - np.roll(np.roll(temp, 1, axis=1), -1, axis=0)[1:-1, 1:-1]
            elif i == 7:
                f += - np.roll(np.roll(temp, 1, axis=1), 1, axis=0)[1:-1, 1:-1]
            elif i == 8:
                f += np.roll(np.roll(temp, -1, axis=1), 1, axis=0)[1:-1, 1:-1]
        return f / 12

    def getNabla_psiy(self):
        f = np.zeros((H, W))
        psi_with_block = copy.deepcopy(self.psi)
        psi_with_block[self.block_mask] = psi_wall
        temp = np.hstack((self.left_wall, np.hstack((psi_with_block, self.right_wall))))
        temp = np.vstack((self.top_bottom_wall, np.vstack((temp, self.top_bottom_wall))))
        for i in range(9):
            if i == 2:
                f += 4 * np.roll(temp, -1, axis=0)[1:-1, 1:-1]
            elif i == 4:
                f += -4 * np.roll(temp, 1, axis=0)[1:-1, 1:-1]
            elif i == 5:
                f += np.roll(np.roll(temp, -1, axis=1), -1, axis=0)[1:-1, 1:-1]
            elif i == 6:
                f += np.roll(np.roll(temp, 1, axis=1), -1, axis=0)[1:-1, 1:-1]
            elif i == 7:
                f += - np.roll(np.roll(temp, 1, axis=1), 1, axis=0)[1:-1, 1:-1]
            elif i == 8:
                f += - np.roll(np.roll(temp, -1, axis=1), 1, axis=0)[1:-1, 1:-1]
        return f / 12

    def getNabla_psi2(self):
        f = np.zeros((H, W))
        psi_with_block = copy.deepcopy(self.psi)
        psi_with_block[self.block_mask] = psi_wall
        temp = np.hstack((self.left_wall, np.hstack((psi_with_block, self.right_wall))))
        temp = np.vstack((self.top_bottom_wall, np.vstack((temp, self.top_bottom_wall))))
        for i in range(9):
            if i == 0:
                f += -20 * temp[1:-1, 1:-1]
            elif i == 1:
                f += 4 * np.roll(temp, -1, axis=0)[1:-1, 1:-1]
            elif i == 2:
                f += 4 * np.roll(temp, -1, axis=1)[1:-1, 1:-1]
            elif i == 3:
                f += 4 * np.roll(temp, 1, axis=1)[1:-1, 1:-1]
            elif i == 4:
                f += 4 * np.roll(temp, 1, axis=0)[1:-1, 1:-1]
            elif i == 5:
                f += np.roll(np.roll(temp, -1, axis=1), -1, axis=0)[1:-1, 1:-1]
            elif i == 6:
                f += np.roll(np.roll(temp, 1, axis=1), -1, axis=0)[1:-1, 1:-1]
            elif i == 7:
                f += np.roll(np.roll(temp, 1, axis=1), 1, axis=0)[1:-1, 1:-1]
            elif i == 8:
                f += np.roll(np.roll(temp, -1, axis=1), 1, axis=0)[1:-1, 1:-1]
        return f / 6

    def getF(self, i):
        f = self.f[i][self.mask] - 1.0 / self.mix_tau * (self.f[i][self.mask] - self.feq[i]) + self.F[i]
        return f

    def getG(self, i):
        g = self.g[i][self.mask] - 1.0 / tau * (self.g[i][self.mask] - self.geq[i])
        return g

    """http://phelafel.technion.ac.il/~drorden/project/ZouHe.pdf"""

    def zou_he_boundary_inlet(self):
        ux = u0
        psi_x = self.getNabla_psix()
        psi_y = self.getNabla_psiy()
        mu = self.getMu_plain()
        rho_inlet = 1 / (1 - ux) * (self.f[0][:, 0] + self.f[2][:, 0] + self.f[4][:, 0] + 2 * (
                self.f[3][:, 0] + self.f[6][:, 0] + self.f[7][:, 0]) - psi_x[:, 0] * mu[:, 0] / 2)
        psi_in = 1.0 - (
                self.g[0][:, 0] + self.g[2][:, 0] + self.g[3][:, 0] + self.g[4][:, 0] + self.g[6][:, 0] + self.g[7][
                                                                                                          :, 0])
        for i in range(9):
            if i == 1:
                self.f[i][1:-1, 0] = self.f[3][1:-1, 0] + 1.5 * ux * rho_inlet[1:-1] - psi_x[1:-1, 0] * mu[1:-1, 0] / 6
                self.g[i][1:-1, 0] = self.w[i] * psi_in[1:-1] / (self.w[1] + self.w[5] + self.w[8])
            if i == 5:
                self.f[i][1:-1, 0] = self.f[7][1:-1, 0] - 0.5 * (
                        self.f[2][1:-1, 0] - self.f[4][1:-1, 0]) + 1.0 / 6.0 * ux * rho_inlet[1:-1] - psi_x[1:-1, 0] * mu[1:-1, 0] / 6 - psi_y[1:-1, 0] * mu[1:-1, 0] / 4
                self.g[i][1:-1, 0] = self.w[i] * psi_in[1:-1] / (self.w[1] + self.w[5] + self.w[8])
            if i == 8:
                self.f[i][1:-1, 0] = self.f[6][1:-1, 0] + 0.5 * (
                        self.f[2][1:-1, 0] - self.f[4][1:-1, 0]) + 1.0 / 6.0 * ux * rho_inlet[1:-1] - psi_x[1:-1, 0] * mu[1:-1, 0] / 6 + psi_y[1:-1, 0] * mu[1:-1, 0] / 4
                self.g[i][1:-1, 0] = self.w[i] * psi_in[1:-1] / (self.w[1] + self.w[5] + self.w[8])

        # left bottom corner node
        self.f[1][0, 0] = self.f[3][0, 0]
        self.g[1][0, 0] = self.g[3][0, 0]
        self.f[2][0, 0] = self.f[4][0, 0]
        self.g[2][0, 0] = self.g[4][0, 0]
        self.f[5][0, 0] = self.f[7][0, 0]
        self.g[5][0, 0] = self.g[7][0, 0]
        self.f[6][0, 0] = 0.5 * (rho_inlet[1] - (self.f[0][0, 0] + self.f[1][0, 0] + self.f[2][0, 0]
                                                 + self.f[3][0, 0] + self.f[4][0, 0] + self.f[5][0, 0] + self.f[7][
                                                     0, 0]))
        self.g[6][0, 0] = self.w[6] * (1.0 - (self.g[0][0, 0] + self.g[1][0, 0] + self.g[2][0, 0]
                                              + self.g[3][0, 0] + self.g[4][0, 0] + self.g[5][0, 0] + self.g[7][
                                                  0, 0])) / (self.w[6] + self.w[8])
        self.f[8][0, 0] = self.f[6][0, 0]
        self.g[8][0, 0] = self.g[6][0, 0]

        # left top corner node
        self.f[1][-1, 0] = self.f[3][-1, 0]
        self.g[1][-1, 0] = self.g[3][-1, 0]
        self.f[4][-1, 0] = self.f[2][-1, 0]
        self.g[4][-1, 0] = self.g[2][-1, 0]
        self.f[8][-1, 0] = self.f[6][-1, 0]
        self.g[8][-1, 0] = self.g[6][-1, 0]
        self.f[5][-1, 0] = 0.5 * (rho_inlet[-2] - (self.f[0][-1, 0] + self.f[1][-1, 0] + self.f[2][-1, 0]
                                                   + self.f[3][-1, 0] + self.f[4][-1, 0] + self.f[6][-1, 0] + self.f[8][
                                                       -1, 0]))
        self.g[5][-1, 0] = self.w[5] * (1.0 - (self.g[0][-1, 0] + self.g[1][-1, 0] + self.g[2][-1, 0]
                                               + self.g[3][-1, 0] + self.g[4][-1, 0] + self.g[6][-1, 0] + self.g[8][
                                                   -1, 0])) / (self.w[5] + self.w[7])
        self.f[7][-1, 0] = self.f[5][-1, 0]
        self.g[7][-1, 0] = self.g[5][-1, 0]

    """Lattice Boltzmann method: fundamentals and engineering applications with computer codes"""
    # open_boundary_with_force
    def zou_he_boundary_outlet(self):
        ux = u0
        psi_x = self.getNabla_psix()[:, -1]
        psi_y = self.getNabla_psiy()[:, -1]
        mu = self.getMu_plain()[:, -1]
        rho_outlet = 1 / (1 + ux) * (self.f[0][:, -1] + self.f[2][:, -1] + self.f[4][:, -1] + 2 * (
                self.f[1][:, -1] + self.f[5][:, -1] + self.f[8][:, -1]) + psi_x * mu / 2)
        psi_out = -1.0 - (self.g[0][:, -1] + self.g[1][:, -1] + self.g[2][:, -1] + self.g[4][:, -1] + self.g[5][:, -1] +
                          self.g[8][:, -1])
        self.f[3][:, -1] = self.f[1][:, -1] - 1.5 * ux * rho_outlet + psi_x * mu / 6
        self.g[3][:, -1] = self.w[3] * psi_out / (self.w[3] + self.w[6] + self.w[7])
        self.f[6][:, -1] = self.f[8][:, -1] - 0.5 * (
                self.f[2][:, -1] - self.f[4][:, -1]) - 1.0 / 6.0 * ux * rho_outlet + psi_y * mu / 4 + psi_x * mu / 6
        self.g[6][:, -1] = self.w[6] * psi_out / (self.w[3] + self.w[6] + self.w[7])
        self.f[7][:, -1] = self.f[5][:, -1] + 0.5 * (
                self.f[2][:, -1] - self.f[4][:, -1]) - 1.0 / 6.0 * ux * rho_outlet - psi_y * mu / 4 + psi_x * mu / 6
        self.g[7][:, -1] = self.w[7] * psi_out / (self.w[3] + self.w[6] + self.w[7])

        self.f[2][0, -1] = self.f[4][0, -1]
        self.g[2][0, -1] = self.g[4][0, -1]
        self.f[4][-1, -1] = self.f[2][-1, -1]
        self.g[4][-1, -1] = self.g[2][-1, -1]


def stream(f, g):
    # with open_boundary on right side
    for i in range(9):
        if i == 1:
            f[i] = np.roll(f[i], 1, axis=1)
            g[i] = np.roll(g[i], 1, axis=1)
        elif i == 2:
            f[i] = np.roll(f[i], 1, axis=0)
            g[i] = np.roll(g[i], 1, axis=0)
        elif i == 3:
            f[i] = np.roll(f[i], -1, axis=1)
            g[i] = np.roll(g[i], -1, axis=1)
        elif i == 4:
            f[i] = np.roll(f[i], -1, axis=0)
            g[i] = np.roll(g[i], -1, axis=0)
        elif i == 5:
            f[i] = np.roll(np.roll(f[i], 1, axis=1), 1, axis=0)
            g[i] = np.roll(np.roll(g[i], 1, axis=1), 1, axis=0)
        elif i == 6:
            f[i] = np.roll(np.roll(f[i], -1, axis=1), 1, axis=0)
            g[i] = np.roll(np.roll(g[i], -1, axis=1), 1, axis=0)
        elif i == 7:
            f[i] = np.roll(np.roll(f[i], -1, axis=1), -1, axis=0)
            g[i] = np.roll(np.roll(g[i], -1, axis=1), -1, axis=0)
        elif i == 8:
            f[i] = np.roll(np.roll(f[i], 1, axis=1), -1, axis=0)
            g[i] = np.roll(np.roll(g[i], 1, axis=1), -1, axis=0)

def update(i, x, y, cc):
    print(i)
    plt.cla()
    plt.pcolor(x, y, cc[i], label='MAX_T{}_Pe{}_M{}_Ca{}_wall{}'.format(MAX_T, Pe, M, Ca, psi_wall), cmap='RdBu')
    # plt.clim(0,1)
    plt.legend()


def bottom_top_wall(f_behind, g_behind, f, g):
    for i in range(9):
        if i == 2:
            f[i][0, :] = f_behind[4][0, :]
            g[i][0, :] = g_behind[4][0, :]
        elif i == 4:
            f[i][-1, :] = f_behind[2][-1, :]
            g[i][-1, :] = g_behind[2][-1, :]
        elif i == 5:
            f[i][0, :] = f_behind[7][0, :]
            g[i][0, :] = g_behind[7][0, :]
        elif i == 6:
            f[i][0, :] = f_behind[8][0, :]
            g[i][0, :] = g_behind[8][0, :]
        elif i == 7:
            f[i][-1, :] = f_behind[5][-1, :]
            g[i][-1, :] = g_behind[5][-1, :]
        elif i == 8:
            f[i][-1, :] = f_behind[6][-1, :]
            g[i][-1, :] = g_behind[6][-1, :]

def main():
    rect_corner_list = []

    # for rectangle block
    # while True:
    #     if count * 20 > 360:
    #         break
    #     if flag:
    #         rect_corner_list.append(((count * 20, 60), ((count + 1) * 20, 80)))
    #         rect_corner_list.append(((count * 20, 140), ((count + 1) * 20, 160)))
    #         rect_corner_list.append(((count * 20, 220), ((count + 1) * 20, 240)))
    #         rect_corner_list.append(((count * 20, 300), ((count + 1) * 20, 320)))
    #         flag = False
    #         count += 2
    #         continue
    #     elif not flag:
    #         rect_corner_list.append(((count * 20, 20), ((count + 1) * 20, 40)))
    #         rect_corner_list.append(((count * 20, 100), ((count + 1) * 20, 120)))
    #         rect_corner_list.append(((count * 20, 180), ((count + 1) * 20, 200)))
    #         rect_corner_list.append(((count * 20, 260), ((count + 1) * 20, 280)))
    #         rect_corner_list.append(((count * 20, 340), ((count + 1) * 20, 360)))
    #         flag = True
    #         count += 2
    #block_psi_all, corner_list = setblock(rect_corner_list)
    cr = Createblock(H, W)
    bb = Bounce_back(H, W)
    circle_list = []
    r = 20

    count = 3
    flag = True
    mabiki = MAX_T // 150

    while True:
        if count * r > 380:
            break
        circle_list.append(((count * r, 2 * r), r))
        circle_list.append(((count * r, 6 * r), r))
        circle_list.append(((count * r, 10 * r), r))
        circle_list.append(((count * r, 14 * r), r))
        circle_list.append(((count * r, 18 * r), r))
        count += 4


    block_psi_all, side_list, concave_list, convex_list = cr.setCirleblock(circle_list)
    block_mask = np.where(block_psi_all == 1, True, False)
    mask = np.logical_not(block_mask)
    cm = Compute(mask)
    cc = np.array([cm.psi])
    for i in range(MAX_T):
        for j in range(9):
            cm.F[j] = cm.getLarge_F(j)
            cm.feq[j] = cm.getfeq(j)
            cm.geq[j] = cm.getgeq(j)
            cm.f[j][mask] = cm.getF(j)
            cm.g[j][mask] = cm.getG(j)
        if i % mabiki == 0:
            cc = np.append(cc, np.array([cm.psi]), axis=0)
            print("timestep:{}".format(i))
        f_behind = copy.deepcopy(cm.f)
        g_behind = copy.deepcopy(cm.g)
        stream(cm.f, cm.g)
        bb.halfway_bounceback_circle(side_list, concave_list, convex_list, f_behind, g_behind, cm.f, cm.g)
        #halfway_bounceback(corner_list, f_behind, g_behind, cm.f, cm.g)
        bottom_top_wall(f_behind[:, 1:-1], g_behind[:, 1:-1], cm.f[:, 1:-1], cm.g[:, 1:-1])
        cm.zou_he_boundary_inlet()
        cm.zou_he_boundary_outlet()
        cm.rho = cm.getRho()
        cm.udpatePsi()
        cm.nabla_psix = cm.getNabla_psix()
        cm.nabla_psiy = cm.getNabla_psiy()
        cm.nabla_psi2 = cm.getNabla_psi2()
        cm.mu = cm.getMu()
        cm.ux = cm.getUx()
        cm.uy = cm.getUy()
        cm.p = cm.getP()
        cm.mix_tau = cm.getMix_tau()
    y = [i for i in range(H)]
    x = [i for i in range(W)]
    fig = plt.figure()
    plt.colorbar(plt.pcolor(x, y, cc[0], cmap='RdBu'))
    ani = animation.FuncAnimation(fig, update, fargs=(x, y, cc), frames=int(len(cc)))
    ani.save('../movies/MAX_T{}_Pe{}_M{}_Ca{}_wall{}.mp4'.format(MAX_T, Pe, M, Ca, psi_wall), fps=10)
    plt.figure()
    plt.pcolor(x, y, cm.psi, label='MAX_T{}_Pe{}_M{}_Ca{}_wall{}'.format(MAX_T, Pe, M, Ca, psi_wall), cmap='RdBu')
    plt.colorbar()
    plt.legend()
    #plt.grid()
    plt.show()
    # plt.savefig('../images/MAX_T{}_Pe{}_M{}_Ca{}_wall{}.png'.format(MAX_T, Pe, M, Ca, psi_wall))


if __name__ == '__main__':
    np.seterr(all='raise')
    try:
        t1 = time.time()
        main()
        t2 = time.time()
        print((t2 - t1) / 60)
    except Warning as e:
        print(e)
```