19
17

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

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

Last updated at Posted at 2018-12-23

#格子ボルツマン法について
研究で使う機会があり, 苦しかったので情報を共有できたらと思います。
格子ボルツマン法(Lattice boltzmann method 以下LBMと略記)はボルツマン方程式を速度, 時間, 座標に関して離散化して解くCFDの一種です。
以下が簡単なシミュレーション方法になります。
Picture2.png

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

おすすめの参考資料

この本が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

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

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)
19
17
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
19
17

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?