Poisson Equation (Steady-State Diffusive Equation)

The Poisson is partial differential equation employed to represent problems as the pressure distribution over steady-state in space. It is characterized by a steady-state version of the difusive equation, i.e, the parameter is independent of time. The partial differential equation is given by

(25)\[ \partial_{\alpha}\left( \nu \partial_{\alpha}p \right) = 0 \]

where \(p\) is the pressure and \(\nu\) diffusive control paremeter.

Numerical Discretization

../_images/mesh-1-diffusive-steady.svg — Fig. 4 1D regular grid example for diffusive equation.

FDM Discretization

Discretizing the Eq. (18) through finite diference method, we have for forward time and central space (FTCS) discretizations:

\[ \partial_{\alpha}\left( \nu \partial_{\alpha}p \right) = \nu\frac{ p(x - \Delta_{x},t) -2p(x,t) + p(x + \Delta_{x},t) }{\Delta_{x}^{2}} =0, \]

for \(\nu\) constant and isolating the term \(u(x,t + \Delta_{t})\), we can describe the \(p\) time evolution:

\[ p(x,t) = \frac{ p(x - \Delta_{x},t) + p(x + \Delta_{x},t) }{2} . \]

Lattice Boltzmann Equation

Describing the problem through the BGK lattice Boltzmann equation (BGK-LB) for the lattice D1Q3:

(26)\[ f_i( x_{\alpha} + c_{i,\alpha} \delta t, t'+\delta t) = f_i(x_{\alpha}, t') -\left( \frac{ f_{i} - f_{i}^{eq} }{ \tau } \right), \]

the equilibrium distribution function is defined by

\[\begin{split} f^{eq}_{i} = \left\{ \begin{array}{ll} \left(w_{0} - 1\right)\widehat{p}(x_{\alpha}, t') & i=0 \\ w_{i}\widehat{p}(x_{\alpha}, t') & i\neq 0 \end{array} \right. , \end{split}\]

where \(\widehat{p}\) is a shifted pressure, shifted by a reference constant pressure \(p_{0}\), i.e., \(p=p_{0}+\widehat{p}\) (This assumption is necessary due to pressure by definition can not be negative, but negative pressure values can be found depending of initial simulation set). The equilibrium moments are given by

(27)\[ \displaystyle\sum_{i=0} f_{i}^{eq}=0, \quad \quad \frac{1}{1-w_{0}}\displaystyle\sum_{i=1} f_{i}^{eq} = \frac{1}{1-w_{0}} \displaystyle\sum_{i=1} f_{i}=\widehat{p}, \quad \quad \displaystyle\sum_{i}f_{i}^{eq} c_{i,\alpha} =0 \quad \quad \textrm{and} \quad \quad \displaystyle\sum_{i}f_{i}^{eq} c_{i,\alpha} c_{i,\beta} =c_{s}^{2} \widehat{p} \delta_{\alpha\beta} . \]

Through Chapman–Enskog analysis, it is demonstrated that the first-order non-equilibrium moment describes the pressure gradient and, consequently, the average velocity:

\[ m^{neq}_{\alpha} = \displaystyle\sum_{i} c_{i,\alpha}\left( f_i-f_i^{eq}\right) = -(c_{s}^{2} \tau \delta_{t})\partial_{\alpha}\widehat{p} \quad \quad \textrm{and} \quad\quad \tau = \displaystyle\frac{\nu}{c_{s}^{2} \delta_{t}} +\frac{1}{2}. \]

Lattice Direction Moments

\[\begin{split}\displaystyle \underbrace{\sum_{i=0} f_{i}^{eq} =\sum_{i=0} f_{i} }_{\textrm{Zero-Order Moment}} =0,\quad \quad \underbrace{\frac{1}{1-w_{0}}\sum_{i=1} f_{i}^{eq} = \frac{1}{1-w_{0}}\sum_{i=1} f_{i} }_{\textrm{Zero-Order Moment}} =\widehat{p},\quad \quad \underbrace{\sum_{i=0} f_{i}^{eq}c_{i,x} }_{\textrm{x-First-Order Moment}} =0,\quad \quad \underbrace{\sum_{i=0} f_{i}^{eq}c_{i,y} }_{\textrm{y-First-Order Moment}} =0,\\ \quad \quad \underbrace{\sum_{i=0} f_{i}^{eq}c_{i,x}c_{i,x} }_{\textrm{xx-Second-Order Moment}} =\frac{\widehat{p}}{3},\quad \quad \underbrace{\sum_{i=0} f_{i}^{eq}c_{i,x}c_{i,y} }_{\textrm{xy-Second-Order Moment}} =0 \quad \quad and \quad \quad \underbrace{\sum_{i=0} f_{i}^{eq}c_{i,y}c_{i,y} }_{\textrm{yy-Second-Order Moment}} =\frac{\widehat{p}}{3}\end{split}\]

Chapmann-Enskog Analysis

Applying the Chapmann-Enskog procedure to LB equation, we expand the term \(f_{i}\left(\boldsymbol{x}+\boldsymbol{e}_i \Delta t, t+\Delta t\right)\) in a Taylor series to recover the derivative form of the equation, i.e.,

(28)\[ f_{i}\left(x_{\alpha} + c_{i,\alpha} \Delta t, t+\Delta t\right)- f_{i}(\boldsymbol{x}, t)=\displaystyle\sum_{j=1}^{\infty}\frac{\Delta t^{j}}{j!}D_{t}^{j}f_{i}= -\left( \frac{ f_{i} - f_{i}^{eq} }{ \tau } \right). \]

Rescaling the dimensionless form of the Eq. (28) in terms of the Knudsen number (\(Kn\)), we have

\[ \displaystyle \sum^{\infty}_{j=1} \frac{Kn^{(j-1)}}{j!} D_{t}^{j} f_{i} = - \frac{1}{Kn}\left( \frac{ f_{i} - f_{eq,i} }{ \tau } \right) \quad \quad \rightarrow \quad \quad \displaystyle \sum^{\infty}_{j=1} \frac{Kn^{(j)}}{j!} D_{t}^{j} f_{i} = - \frac{ f_{i} - f_{i}^{eq} }{ \tau }, \]

applying the asymptotic expansion in both the distribution function (\(f_{i}=f_{i}^{(0)}+Kn f_{i}^{(1)} + Kn^{2} f_{i}^{(2)}+\cdots\)) and time partial derivative (\(\partial_{t}=\partial_{t}^{(0)}+ Kn \partial_{t}^{(1)}+Kn^{2} \partial_{t}^{(2)}+\cdots\)), and separating the equation in orders up to the order \(Kn^{2}\):

(29)\[\begin{split} \begin{array}{ll} (Kn^{(0)}):& f_{i}^{(0)}= f_{i}^{eq},\\ (Kn^{(1)}):& \left( \partial_{t}^{(0)} + c_{i,\alpha}\partial_{\alpha} \right)f_{i}^{(0)}= - \displaystyle\frac{ f_{i}^{(1)} }{ \tau } , \\ (Kn^{(2)}):& \partial_{t}^{(1)} f_{i}^{(0)} + \left( \partial_{t}^{(0)} + c_{i,\alpha}\partial_{\alpha} \right)f_{i}^{(0)} + \displaystyle\frac{\left( \partial_{t}^{(0)} + c_{i,\alpha}\partial_{\alpha} \right)^{2} f_{i}^{(0)}}{2}= - \displaystyle\frac{ f_{i}^{(2)} }{ \tau },\\ \textrm{ou}\\ (Kn^{(2)}):& \partial_{t}^{(1)} f_{i}^{(0)} + \left( \partial_{t}^{(0)} + c_{i,\alpha}\partial_{\alpha} \right)\displaystyle\left(1 - \frac{1}{2\tau}\right) f_{i}^{(1)} = - \displaystyle\frac{ f_{i}^{(2)} }{ \tau } \quad\quad \rightarrow \quad\quad f_{i}^{(1)}\textrm{ Formulation} \\ \textrm{ou}\\ (Kn^{(2)}):& \partial_{t}^{(1)} f_{i}^{(0)} + \left( \partial_{t}^{(0)} + c_{i,\alpha}\partial_{\alpha} \right)^{2} \displaystyle\left(\frac{1}{2} - \tau\right) f_{i}^{(0)} = - \displaystyle\frac{ f_{i}^{(2)} }{ \tau } \quad\quad \rightarrow \quad\quad f_{i}^{(0)}\textrm{ Formulation}. \end{array} \end{split}\]

Zero-Order Moment Balance

To retrieve the balance equation, we sum the Eq. (29) for \(Kn^{(1)}\) and \(Kn^{(2)}\) over \(\sum_{i=0} \):

\[\begin{split} \begin{array}{l} (Kn^{(1)}): \displaystyle\sum_{i}\left( \partial_{t}^{(0)} + c_{i,\alpha}\partial_{\alpha} \right)f_{i}^{(0)}= \displaystyle\sum_{i} \left(- \displaystyle\frac{ f_{i}^{(1)} }{ \tau } \right),\\ (Kn^{(1)}): 0 =0, \end{array} \end{split}\]

and

\[\begin{split} \begin{array}{l} (Kn^{(2)}): \displaystyle\sum_{i}\left( \partial_{t}^{(1)} f_{i}^{(0)} + \left( \partial_{t}^{(0)} + c_{i,\alpha}\partial_{\alpha} \right)^{2} \displaystyle\left( \frac{1}{2}-\tau\right) f_{i}^{(0)} \right) = \displaystyle\sum_{i}\left( - \displaystyle\frac{ f_{i}^{(2)} }{ \tau } \right), \\ (Kn^{(2)}): \left( \displaystyle\frac{1}{2}-\tau\right) \partial_{\alpha}\partial_{\beta} \left(c_{s}^{2} \widehat{p} \delta_{\alpha\beta} \right) =0. \end{array} \end{split}\]

or

(30)\[\begin{split} \begin{array}{l} (Kn^{(2)}):& \displaystyle\sum_{i}\left[ \partial_{t}^{(1)} f_{i}^{(0)} - \left(1-\displaystyle\frac{1}{2\tau}\right)\left( \partial_{t}^{(0)} + c_{i,\alpha}\partial_{\alpha} \right) f_{i}^{(1)} \right]= \displaystyle\sum_{i}\left( - \displaystyle\frac{f_{i}^{(2)}}{\tau} \right), \nonumber \\ (Kn^{(2)}):& - \left(1-\displaystyle\frac{1}{2\tau}\right) m^{(1)}_{\alpha} = 0. \end{array} \end{split}\]

where \(m^{(1)}_{\alpha}\) is the first-order moment of \(f_{i}^{(1)}\) and \(1/(1-w_{0})\sum_{i=1}f_{i}^{(j)}=0\) for \(j\geq 1\) due to imposition of \(\widehat{p}\) conservation. To compute the moment \(m^{(1)}_{\alpha}\), we multiply Eq. (29) by \(c_{i,\alpha}\) and sum over all lattice directions:

(31)\[\begin{split} \begin{array}{l} (Kn^{(1)}): \displaystyle\sum_{i}c_{i,\alpha}\left( \partial_{t}^{(0)} + e_{\beta,i}\partial_{\beta} \right)f_{i}^{(0)}= \displaystyle\sum_{i} c_{i,\alpha} \left( -\frac{f_{i}^{(1)}}{\tau} \right),\\ (Kn^{(1)}): m^{(1)}_{\alpha} =-\tau \partial_{\beta}\left( c_{s}^{2} \widehat{p} \delta_{\alpha\beta} \right). \end{array} \end{split}\]

By substituting Eq. (31) into Eq. (30), we recover the Eq. (25), where \(\nu=c_{s}^{2}(\tau-1/2)\). In the regularized formulation of the lattice Boltzmann equation, the first-order correction is approximated as \(f_{i}^{(1)}\approx (f_{i}-f_{i}^{(0)})=f_{i}^{neq}\), where higher-order Knudsen moments are filtered out and the collision term is reformulated in terms of \(f_{i}^{neq}\).

Boundary Conditions

The boundary conditions for the lattice D2Q5 can be derived by solving a linear system of known moments, as illustrated in Boundary Conditions: Deduction for the Lattice D2Q5 - Poisson Equation. For the D2Q5 lattice, the unknown distribution functions at each boundary are determined as summarized in Table below.

Boundaries	Layers
	North	South
Unknown \(f_{i}\)	\(f_2=\widehat{p}_{n}(1-w_{0}) - (f_1+f_3+f_4)\)	\(f_4=\widehat{p}_{s}(1-w_{0}) - (f_1+f_2+f_3)\)
	East	West
	\(f_3=\widehat{p}_{e}(1-w_{0}) - (f_1+f_2+f_4)\)	\(f_1=\widehat{p}_{w}(1-w_{0}) - (f_2+f_3+f_4)\)
Boundaries	Concave Corners
	Northwest	Northeast
Unknown \(f_{i}\)	\(f_1=-(\partial_{x}\widehat{p}+\partial_{y}\widehat{p})\tau c_{s}^{2}/2 - f_2 + \widehat{p}_{nw}/3\)	\(f_3=(\partial_{x}\widehat{p}-\partial_{y}\widehat{p})\tau c_{s}^{2}/2 - f_2 + \widehat{p}_{ne}/3\)
	\(f_4=(\partial_{x}\widehat{p}+\partial_{y}\widehat{p})\tau c_{s}^{2}/2 - f_3 + \widehat{p}_{nw}/3\)	\(f_4=-(\partial_{x}\widehat{p}-\partial_{y}\widehat{p})\tau c_{s}^{2}/2 - f_1 + \widehat{p}_{ne}/3\)
	Southwest	Southeast
	\(f_1=-(\partial_{x}\widehat{p}-\partial_{y}\widehat{p})\tau c_{s}^{2}/2 - f_4 + \widehat{p}_{sw}/3\)	\(f_2=-(\partial_{x}\widehat{p}+\partial_{y}\widehat{p})\tau c_{s}^{2}/2 - f_1 + \widehat{p}_{se}/3\)
	\(f_2=(\partial_{x}\widehat{p}-\partial_{y}\widehat{p})\tau c_{s}^{2}/2 - f_3 + \widehat{p}_{sw}/3\)	\(f_3=(\partial_{x}\widehat{p}+\partial_{y}\widehat{p})\tau c_{s}^{2}/2 - f_4 + \widehat{p}_{se}/3\)

1st Benchmark: Single sine mode on a finite rod (Dirichlet)

PDE: \((\phi_t = \partial_{x}(\nu\partial_{x}\phi), (0<x<L,\ t>0)\)

BC: \((\phi(0,t)=\phi(L,t)=0)\)

IC: \(\phi(x,0)=\sin \left(\frac{\pi x}{L}\right)\)

Solution

\[ \phi(x,t)=\sin\left(\frac{\pi x}{L}\right) \exp \left(-\nu \frac{\pi^2}{L^2} t\right). \]

FDM Solution:

LBM D1Q3 Solution With Steady-Scheme:

Show code cell source Hide code cell source

import warnings
warnings.filterwarnings("ignore")
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams['mathtext.fontset'] = 'cm'
#--------------------------------------- Parameters (physical units) ------------------------------------------
L = 1.0                # Length of the domain
nu = 0.1/3.0           # Diffusion coefficient
t_end = 2.0                # Total simulation time
#--------------------------------------- Lattice-Properties-D1Q3 ----------------------------------------------
cs=1.0/np.sqrt(3.0);
w = np.array([4.0/6.0, 1.0/6.0, 1.0/6.0],dtype="float64")
cx = np.array([0, 1, -1],dtype="int8")  
#------------------------------------- Parameters (numerical units) -------------------------------------------
Nx=21                   #Square Domain Length
dx = L / (Nx-1)         # Spatial step size
c=1.0*2**(2)            # c=dx/dt
dt=dx/c                 # Time step size
nt = int(t_end / dt)        # Time step number
nue=nu/(c*dx)       # Diffusion coefficient
tau=nue/cs**2+0.5   # Relaxation Time
print(f"dx={dx:.4e}\t dt={dt:.4e}\t nt={nt:d}") # Print values for check
print(f"nue={nue:.4f}\t tau={tau:.4f}")
#--------------------------------- Initialization - Numerical Arrays -----------------------------------------
T = np.zeros((Nx),dtype="float64")
x = np.linspace(0, Nx-1, Nx)
T = np.sin(np.pi * x / (Nx-1))
f = np.zeros((3,Nx),dtype="float64")
fp = np.zeros((3,Nx),dtype="float64")
AvT = np.zeros((nt),dtype="float64")
for k in range(0,3):
    fp[k,:]=w[k]*T[:]
#----------------------------------------- Maind Loop --------------------------------------------------------
for t in range(nt):
    #-----------------streaming-------------------
    for k in range(0,3):
        f[k,:]=np.roll(fp[k,:], cx[k], axis=0)
    #-----------------Boundaries-----------------------
    f[1,0]= 0.0 - f[0,0]-f[2,0]
    f[2,Nx-1]= 0.0 - f[0,Nx-1]-f[1,Nx-1]
    #----------------------Macro------------------
    T[:] = (f[1,:]+f[2,:])/(1-w[0])
    AvT[t]=np.sum(T)/Nx;
    #--------------------Collision----------------
    fp[0,:]=f[0,:] - (f[0,:] - (w[0]-1)*T[:])/tau
    fp[1,:]=f[1,:] - (f[1,:] - w[1]*T[:])/tau
    fp[2,:]=f[2,:] - (f[2,:] - w[2]*T[:])/tau

LBM D1Q2 Solution (Not Accept the Steady-Scheme):

Show code cell source Hide code cell source

import warnings
warnings.filterwarnings("ignore")
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams['mathtext.fontset'] = 'cm'
#--------------------------------------- Parameters (physical units) ------------------------------------------
L = 1.0                # Length of the domain
nu = 0.1/2.0           # Diffusion coefficient
t_end = 2.0            # Total simulation time
#--------------------------------------- Lattice-Properties-D1Q2 ----------------------------------------------
cs=1.0/np.sqrt(2.0);
w = np.array([1.0/2.0, 1.0/2.0],dtype="float64")
cx = np.array([1, -1],dtype="int8")  
#------------------------------------- Parameters (numerical units) -------------------------------------------
Nx=21                   # Square Domain Length
dx = L / (Nx-1)         # Spatial step size
c=1.0*2**(2)            # c=dx/dt
dt=dx/c                 # Time step size
nt = int(t_end / dt)        # Time step number
nue=nu/(c*dx)           # Diffusion coefficient
tau=nue/cs**2+0.5       # Relaxation Time
print(f"dx={dx:.4e}\t dt={dt:.4e}\t nt={nt:d}") # Print values for check
print(f"nue={nue:.4f}\t tau={tau:.4f}")
#--------------------------------- Initialization - Numerical Arrays -----------------------------------------
T2 = np.zeros((Nx),dtype="float64")
x = np.linspace(0, Nx-1, Nx)
T2 = np.sin(np.pi * x / (Nx-1))
f = np.zeros((2,Nx),dtype="float64")
fp = np.zeros((2,Nx),dtype="float64")
AvT2 = np.zeros((nt),dtype="float64")
fp[0,:]=w[0]*T2[:]
fp[1,:]=w[1]*T2[:]
#----------------------------------------- Maind Loop --------------------------------------------------------
for t in range(nt):
    #-----------------streaming-------------------
    f[0,:]=np.roll(fp[0,:], cx[0], axis=0)
    f[1,:]=np.roll(fp[1,:], cx[1], axis=0)
    #-----------------Boundaries-----------------------
    f[0,0]= 0.0 - f[1,0]
    f[1,Nx-1]= 0.0 - f[0,Nx-1]
    #----------------------Macro------------------
    T2[:] = f[0,:]+f[1,:]
    AvT2[t]=np.sum(T2)/Nx;
    #--------------------Collision----------------
    fp[0,:]=f[0,:] - (f[0,:] - w[0]*T2[:])/tau
    fp[1,:]=f[1,:] - (f[1,:] - w[1]*T2[:])/tau

LBM D1Q3 Solution (Transient):

Show code cell source Hide code cell source

import warnings
warnings.filterwarnings("ignore")
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams['mathtext.fontset'] = 'cm'
#--------------------------------------- Parameters (physical units) ------------------------------------------
L = 1.0                # Length of the domain
nu = 0.1/3.0           # Diffusion coefficient
t = 2.0                # Total simulation time
#--------------------------------------- Lattice-Properties-D1Q3 ----------------------------------------------
cs=1.0/np.sqrt(3.0);
w = np.array([4.0/6.0, 1.0/6.0, 1.0/6.0],dtype="float64")
cx = np.array([0, 1, -1],dtype="int8")
#------------------------------------- Parameters (numerical units) -------------------------------------------
Nx=21                   # Square Domain Length
dx = L / (Nx-1)         # Spatial step size
c=1.0*2**(2)            # c=dx/dt
dt=dx/c                 # Time step size
nt = int(t / dt)        # Time step number
nue=nu/(c*dx)       # Diffusion coefficient
tau=nue/cs**2+0.5   # Relaxation Time
print(f"dx={dx:.4e}\t dt={dt:.4e}\t nt={nt:d}") # Print values for check
print(f"nue={nue:.4f}\t tau={tau:.4f}")
#--------------------------------- Initialization - Numerical Arrays -----------------------------------------
Tt = np.zeros((Nx),dtype="float64")
x = np.linspace(0, Nx-1, Nx)
Tt = np.sin(np.pi * x / (Nx-1))
f = np.zeros((3,Nx),dtype="float64")
fp = np.zeros((3,Nx),dtype="float64")
AvTt = np.zeros((nt),dtype="float64")
for k in range(0,3):
    fp[k,:]=w[k]*Tt[:]
#----------------------------------------- Maind Loop --------------------------------------------------------
for t in range(nt):
    #-----------------streaming-------------------
    for k in range(0,3):
        f[k,:]=np.roll(fp[k,:], cx[k], axis=0)
    #-----------------Boundaries-----------------------
    f[1,0]= 0.0 - f[0,0]-f[2,0]
    f[2,Nx-1]= 0.0 - f[0,Nx-1]-f[1,Nx-1]
    #----------------------Macro------------------
    Tt[:]=f[0,:]+f[1,:]+f[2,:]
    AvTt[t]=np.sum(Tt)/Nx;
    #--------------------Collision----------------
    for k in range(0,3):
        fp[k,:]=f[k,:] - (f[k,:] - w[k]*Tt[:])/tau

Numerical and Analytical Comparisson:

../_images/5373282b09696c9b84b7ab476e231d8c643eadbd5a247eec673bcb11d40e233a.png

2nd Benchmark: 2D Poisson Pressure Field

Analyzing the accuracy of Poisson equation solver, we apply it to a two-dimensional analytical case. The geometry is described by a square domain of length \(L\) initialized with a constant dimensionless pressure \(\widehat{p}(x,y,t=0)=0\). The boundary conditions of the geometry are given by

\[ \widehat{p}(x=0,y,t)=\displaystyle\frac{sinh\left(\pi(1-y)\right)}{sinh\left(\pi\right)}, \quad \quad \widehat{p}(x=L,y,t)= -\frac{sinh\left(\pi(1-y)\right)}{sinh\left(\pi\right)}, \quad \quad \widehat{p}(x,y=0,t)=cos(\pi x), \quad \quad \widehat{p}(x,y=L,t)=0. \]

The pressure analytical solution for Eq. (25) pressure is given

\[ \widehat{p}(x,y,t)=cos(\pi x) \displaystyle\frac{sinh\left(\pi(1-y)\right)}{sinh\left(\pi\right)} \]

where \(Dp_{0}\) is considered equal to 1, consequently, velocity values can be analytically obtained by the equation:

\[ u_{x}=D\partial_{x}p= \pi sin(\pi x) \displaystyle\frac{sinh\left(\pi(1-y)\right)}{sinh\left(\pi\right)}, \quad \quad u_{y}=D\partial_{y}p= \pi cos(\pi x) \displaystyle\frac{cosh\left(\pi(1-y)\right)}{sinh\left(\pi\right)}. \]

LBM D2Q5 Code

Show code cell source Hide code cell source

import time
start=time.time()
for kk in range(0,mstep):
    #*******Streaming********
    for k in range(0,5):
        f[k,:,:]=np.roll(np.roll(fp[k,:,:], cx[k], axis=0), cy[k], axis=1)
    #************Boundary-Conditions**********
    f[4,1:-1,Ny-1]=-(f[1,1:-1,Ny-1]+f[2,1:-1,Ny-1]+f[3,1:-1,Ny-1]) #North
    f[2,1:-1,0]=np.cos(np.pi*x[1:-1])*(1.0-w[0])-(f[1,1:-1,0]+f[3,1:-1,0]+f[4,1:-1,0]) #South
    f[1,0,1:-1]=np.sinh(np.pi*(1-y[1:-1]))/np.sinh(np.pi)*(1.0-w[0])-(f[2,0,1:-1]+f[3,0,1:-1]+f[4,0,1:-1]) #West
    f[3,Nx-1,1:-1]=-np.sinh(np.pi*(1-y[1:-1]))/np.sinh(np.pi)*(1.0-w[0])-(f[1,Nx-1,1:-1]+f[2,Nx-1,1:-1]+f[4,Nx-1,1:-1]) #East
    #************Corner-Boundaries**********
    #----------------NorthWest----------------------
    f[1,0,Ny-1]=-(dxpn[0]+dypw[Ny-1])*tau*cs**2/2 -f[2,0,Ny-1]
    f[4,0,Ny-1]=(dxpn[0]+dypw[Ny-1])*tau*cs**2/2 -f[3,0,Ny-1]
    #----------------NorthEast----------------------
    f[3,Nx-1,Ny-1]=(dxpn[Nx-1]-dype[Ny-1])*tau*cs**2/2 -f[2,Nx-1,Ny-1]
    f[4,Nx-1,Ny-1]=-(dxpn[Nx-1]-dype[Ny-1])*tau*cs**2/2 -f[1,Nx-1,Ny-1]
    #----------------SouthWest----------------------
    f[1,0,0]=-(dxps[0]-dypw[0])*tau*cs**2/2 -f[4,0,0] +np.cos(np.pi*x[0])/3.0
    f[2,0,0]=(dxps[0]-dypw[0])*tau*cs**2/2 -f[3,0,0] +np.cos(np.pi*x[0])/3.0
    #----------------SouthEast----------------------
    f[2,Nx-1,0]=-(dxps[Nx-1]+dype[0])*tau*cs**2/2 -f[1,Nx-1,0] +np.cos(np.pi*x[Nx-1])/3.0
    f[3,Nx-1,0]=(dxps[Nx-1]+dype[0])*tau*cs**2/2 -f[4,Nx-1,0] +np.cos(np.pi*x[Nx-1])/3.0
    #***********Macro********
    p=(f[1,:,:]+f[2,:,:]+f[3,:,:]+f[4,:,:])/(1.0-w[0])
    #*******Collision********
    fp[0,:,:]=f[0,:,:]-(f[0,:,:]-(w[0]-1.0)*p[:,:])/tau # Collision Lattice Index 0
    for k in range(1,5):
        fp[k,:,:]=f[k,:,:]-(f[k,:,:]-w[k]*p[:,:])/tau # Collision Lattice Index 1-4

end=time.time();
runtime = end-start;
nodes_updated = mstep*Nx*Ny;
speed = nodes_updated/(1e6*runtime);
print('runtime=',runtime)
print('Mlups=',speed)

../_images/51eb6f083818aba1daf6bf87e7a2372c3166f164096041a7961a27338f443aac.png

../_images/79a3026090fb6832c37ddd37107862b1157838ae39be8cbcf4ad05ba9f3ec2d4.png

../_images/e949468fcf70e72ed79b084de71fc241ef73eadd8565c33e36c578d549702556.png

../_images/c3190a2e88b0477185df0eba7d063926eafd0d44de499205ddf6a63507b0e244.png