Single_phase.lbm_solver_3d
This file is the non-objective oriented version of singlephase solver without using class. At the begining of the this file it define some variable first.
#import some package
import taichi as ti
import numpy as np
from pyevtk.hl import gridToVTK
import time
#initialize taichi with cpu, dunamic index, disable profiler and disables printing the intermediate representation
ti.init(arch=ti.cpu, dynamic_index=True, kernel_profiler=False, print_ir=False)
#enable projection
enable_projection = True
#nx,ny,nz = 100,50,5
#define 131x131x131 and zero external force
nx,ny,nz = 131,131,131
fx,fy,fz = 0.0e-6,0.0,0.0
#viscosity=0.1
niu = 0.1
#Boundary condition mode: 0=periodic, 1= fix pressure, 2=fix velocity; boundary pressure value (rho); boundary velocity value for vx,vy,vz
bc_x_left, rho_bcxl, vx_bcxl, vy_bcxl, vz_bcxl = 1, 1.0, 0.0e-5, 0.0, 0.0 #Boundary x-axis left side
bc_x_right, rho_bcxr, vx_bcxr, vy_bcxr, vz_bcxr = 1, 0.995, 0.0, 0.0, 0.0 #Boundary x-axis right side
bc_y_left, rho_bcyl, vx_bcyl, vy_bcyl, vz_bcyl = 0, 1.0, 0.0, 0.0, 0.0 #Boundary y-axis left side
bc_y_right, rho_bcyr, vx_bcyr, vy_bcyr, vz_bcyr = 0, 1.0, 0.0, 0.0, 0.0 #Boundary y-axis right side
bc_z_left, rho_bczl, vx_bczl, vy_bczl, vz_bczl = 0, 1.0, 0.0, 0.0, 0.0 #Boundary z-axis left side
bc_z_right, rho_bczr, vx_bczr, vy_bczr, vz_bczr = 0, 1.0, 0.0, 0.0, 0.0 #Boundary z-axis right side
#define old density distribution funciton nx*ny*nz*19
f = ti.field(ti.f32,shape=(nx,ny,nz,19))
#define new density distribution function nx*ny*nz*19
F = ti.field(ti.f32,shape=(nx,ny,nz,19))
#define density nx*ny*nz
rho = ti.field(ti.f32, shape=(nx,ny,nz))
#define velocity nx*ny*nz
v = ti.Vector.field(3,ti.f32, shape=(nx,ny,nz))
#define lattice speed 3*19
e = ti.Vector.field(3,ti.i32, shape=(19))
#define s diagonal 19 dimension vector
S_dig = ti.field(ti.f32,shape=(19))
#define another lattice speed 3*19
e_f = ti.Vector.field(3,ti.f32, shape=(19))
#define weight parameter 19 dimesnion vector
w = ti.field(ti.f32, shape=(19))
#define solid flag nx*ny*nz
solid = ti.field(ti.i32,shape=(nx,ny,nz))
#define vector for streaming 19 dimensional vector
LR = ti.field(ti.i32,shape=(19))
#define external force with a 3 dimensional vector
ext_f = ti.Vector.field(3,ti.f32,shape=())
#define velocity in x,y,z direction with 3 dimensional vector
bc_vel_x_left = ti.Vector.field(3,ti.f32, shape=())
bc_vel_x_right = ti.Vector.field(3,ti.f32, shape=())
bc_vel_y_left = ti.Vector.field(3,ti.f32, shape=())
bc_vel_y_right = ti.Vector.field(3,ti.f32, shape=())
bc_vel_z_left = ti.Vector.field(3,ti.f32, shape=())
bc_vel_z_right = ti.Vector.field(3,ti.f32, shape=())
#define transforming matrix 19*19
M = ti.field(ti.f32, shape=(19,19))
#define inverse of transforming matrix
inv_M = ti.field(ti.f32, shape=(19,19))
#define single relaxation parameter
tau_f=3.0*niu+0.5
#define single relaxation frequency
s_v=1.0/tau_f
#define other parameter in the s diagonal
s_other=8.0*(2.0-s_v)/(8.0-s_v)
#define s matrix but not used
S_np = np.zeros((19,19))
S_np[0,0]=0; S_np[1,1]=s_v; S_np[2,2]=s_v; S_np[3,3]=0; S_np[4,4]=s_other; S_np[5,5]=0;
S_np[6,6]=s_other; S_np[7,7]=0; S_np[8,8]=s_other; S_np[9,9]=s_v; S_np[10,10]=s_v; S_np[11,11]=s_v;
S_np[12,12]=s_v; S_np[13,13]=s_v; S_np[14,14]=s_v; S_np[15,15]=s_v; S_np[16,16]=s_other; S_np[17,17]=s_other;
S_np[18,18]=s_other
#define numpy array version of s diagonal.
S_dig_np = np.array([0,s_v,s_v,0,s_other,0,s_other,0,s_other, s_v, s_v,s_v,s_v,s_v,s_v,s_v,s_other,s_other,s_other])
#define numpy version of transforming matrix
M_np = np.array([[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1],
[-1,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1],
[1,-2,-2,-2,-2,-2,-2,1,1,1,1,1,1,1,1,1,1,1,1],
[0,1,-1,0,0,0,0,1,-1,1,-1,1,-1,1,-1,0,0,0,0],
[0,-2,2,0,0,0,0,1,-1,1,-1,1,-1,1,-1,0,0,0,0],
[0,0,0,1,-1,0,0,1,-1,-1,1,0,0,0,0,1,-1,1,-1],
[0,0,0,-2,2,0,0,1,-1,-1,1,0,0,0,0,1,-1,1,-1],
[0,0,0,0,0,1,-1,0,0,0,0,1,-1,-1,1,1,-1,-1,1],
[0,0,0,0,0,-2,2,0,0,0,0,1,-1,-1,1,1,-1,-1,1],
[0,2,2,-1,-1,-1,-1,1,1,1,1,1,1,1,1,-2,-2,-2,-2],
[0,-2,-2,1,1,1,1,1,1,1,1,1,1,1,1,-2,-2,-2,-2],
[0,0,0,1,1,-1,-1,1,1,1,1,-1,-1,-1,-1,0,0,0,0],
[0,0,0,-1,-1,1,1,1,1,1,1,-1,-1,-1,-1,0,0,0,0],
[0,0,0,0,0,0,0,1,1,-1,-1,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,-1,-1],
[0,0,0,0,0,0,0,0,0,0,0,1,1,-1,-1,0,0,0,0],
[0,0,0,0,0,0,0,1,-1,1,-1,-1,1,-1,1,0,0,0,0],
[0,0,0,0,0,0,0,-1,1,1,-1,0,0,0,0,1,-1,1,-1],
[0,0,0,0,0,0,0,0,0,0,0,1,-1,-1,1,-1,1,1,-1]])
#define inverse of transforming matrix using inv function in linalg package
inv_M_np = np.linalg.inv(M_np)
#define index for streaming
LR_np = np.array([0,2,1,4,3,6,5,8,7,10,9,12,11,14,13,16,15,18,17])
#assign numpy version to M.np to M
M.from_numpy(M_np)
#assign numpy version of inverser matrix inv_M_np to inv_M
inv_M.from_numpy(inv_M_np)
#assign numpy versio of LR array to LR
LR.from_numpy(LR_np)
#assign fx,fy,fz to vector external force
ext_f[None] = ti.Vector([fx,fy,fz])
#assign numpy version of S diagnal S_dig_np to S_dig
S_dig.from_numpy(S_dig_np)
#make inv_M,M,LR,S_dig not modified
ti.static(inv_M)
ti.static(M)
ti.static(LR)
ti.static(S_dig)
#create mesh nx*ny*nz
x = np.linspace(0, nx, nx)
y = np.linspace(0, ny, ny)
z = np.linspace(0, nz, nz)
#numpy meshgrid from x,y,z 1d array to 3d array X,Y,Z here use ij indexing
X, Y, Z = np.meshgrid(x, y, z, indexing='ij')
feq(k,rho_local,u) calculate the equilibrium density distribution function in velocity space
# taichi funciton
@ti.func
def feq(k,rho_local, u):
eu = e[k].dot(u)
uv = u.dot(u)
#calculate the equilibrium density distribution function
feqout = w[k]*rho_local*(1.0+3.0*eu+4.5*eu*eu-1.5*uv)
#print(k, rho_local, w[k])
return feqout
init() initialize velocity=0, density=1 and denisty distribution function= equilibrium density distribution function
@ti.kernel
def init():
for i,j,k in rho:
rho[i,j,k] = 1.0
v[i,j,k] = ti.Vector([0,0,0])
for s in range(19):
f[i,j,k,s] = feq(s,1.0,v[i,j,k])
F[i,j,k,s] = feq(s,1.0,v[i,j,k])
#print(F[i,j,k,s], feq(s,1.0,v[i,j,k]))
init_geo() load geometry file
def init_geo(filename):
#load data
in_dat = np.loadtxt(filename)
#reshape it with column major
in_dat = np.reshape(in_dat, (nx,ny,nz),order='F')
return in_dat
static_init() initialize lattixe speed weight parameter and boundary velocity
@ti.kernel
def static_init():
if ti.static(enable_projection): # No runtime overhead
#initialize lattice speed
e[0] = ti.Vector([0,0,0])
e[1] = ti.Vector([1,0,0]); e[2] = ti.Vector([-1,0,0]); e[3] = ti.Vector([0,1,0]); e[4] = ti.Vector([0,-1,0]);e[5] = ti.Vector([0,0,1]); e[6] = ti.Vector([0,0,-1])
e[7] = ti.Vector([1,1,0]); e[8] = ti.Vector([-1,-1,0]); e[9] = ti.Vector([1,-1,0]); e[10] = ti.Vector([-1,1,0])
e[11] = ti.Vector([1,0,1]); e[12] = ti.Vector([-1,0,-1]); e[13] = ti.Vector([1,0,-1]); e[14] = ti.Vector([-1,0,1])
e[15] = ti.Vector([0,1,1]); e[16] = ti.Vector([0,-1,-1]); e[17] = ti.Vector([0,1,-1]); e[18] = ti.Vector([0,-1,1])
#initialize lattice speed
e_f[0] = ti.Vector([0,0,0])
e_f[1] = ti.Vector([1,0,0]); e_f[2] = ti.Vector([-1,0,0]); e_f[3] = ti.Vector([0,1,0]); e_f[4] = ti.Vector([0,-1,0]);e_f[5] = ti.Vector([0,0,1]); e_f[6] = ti.Vector([0,0,-1])
e_f[7] = ti.Vector([1,1,0]); e_f[8] = ti.Vector([-1,-1,0]); e_f[9] = ti.Vector([1,-1,0]); e_f[10] = ti.Vector([-1,1,0])
e_f[11] = ti.Vector([1,0,1]); e_f[12] = ti.Vector([-1,0,-1]); e_f[13] = ti.Vector([1,0,-1]); e_f[14] = ti.Vector([-1,0,1])
e_f[15] = ti.Vector([0,1,1]); e_f[16] = ti.Vector([0,-1,-1]); e_f[17] = ti.Vector([0,1,-1]); e_f[18] = ti.Vector([0,-1,1])
#intialize weight parameter
w[0] = 1.0/3.0; w[1] = 1.0/18.0; w[2] = 1.0/18.0; w[3] = 1.0/18.0; w[4] = 1.0/18.0; w[5] = 1.0/18.0; w[6] = 1.0/18.0;
w[7] = 1.0/36.0; w[8] = 1.0/36.0; w[9] = 1.0/36.0; w[10] = 1.0/36.0; w[11] = 1.0/36.0; w[12] = 1.0/36.0;
w[13] = 1.0/36.0; w[14] = 1.0/36.0; w[15] = 1.0/36.0; w[16] = 1.0/36.0; w[17] = 1.0/36.0; w[18] = 1.0/36.0;
#intialize boundary velocity
bc_vel_x_left[None] = ti.Vector([vx_bcxl, vy_bcxl, vz_bcxl])
bc_vel_x_right[None] = ti.Vector([vx_bcxr, vy_bcxr, vz_bcxr])
bc_vel_y_left[None] = ti.Vector([vx_bcyl, vy_bcyl, vz_bcyl])
bc_vel_y_right[None] = ti.Vector([vx_bcyr, vy_bcyr, vz_bcyr])
bc_vel_z_left[None] = ti.Vector([vx_bczl, vy_bczl, vz_bczl])
bc_vel_z_right[None] = ti.Vector([vx_bczr, vy_bczr, vz_bczr])
multiply_M calculate denisty distribution function in momentum space M*f=m
@ti.func
def multiply_M(i,j,k):
out = ti.Vector([0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0])
for index in range(19):
for s in range(19):
#calculte m=M*f here
out[index] += M[index,s]*F[i,j,k,s]
#print(i,j,k, index, s, out[index], M[index,s], F[i,j,k,s])
return out
GuoF(i,j,k,s,u) calculate Guo’s Force scheme
@ti.func
def GuoF(i,j,k,s,u):
out=0.0
for l in range(19):
#calculate Guo's force here
out += w[l]*((e_f[l]-u).dot(ext_f[None])+(e_f[l].dot(u)*(e_f[l].dot(ext_f[None]))))*M[s,l]
return out
meq_vec(rho_local,u) calculate equilibrium density distribution function in momentum space
@ti.func
def meq_vec(rho_local,u):
out = ti.Vector([0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0])
out[0] = rho_local; out[3] = u[0]; out[5] = u[1]; out[7] = u[2];
out[1] = u.dot(u); out[9] = 2*u.x*u.x-u.y*u.y-u.z*u.z; out[11] = u.y*u.y-u.z*u.z
out[13] = u.x*u.y; out[14] = u.y*u.z; out[15] = u.x*u.z
return out
collison() define the prcoess of collision
@ti.kernel
def colission():
for i,j,k in rho:
#if it is fluid
if (solid[i,j,k] == 0):
#calculate m
m_temp = multiply_M(i,j,k)
#calculate meq
meq = meq_vec(rho[i,j,k],v[i,j,k])
for s in range(19):
#calculate -s*(m-meq)
m_temp[s] -= S_dig[s]*(m_temp[s]-meq[s])
#add Guo's force
m_temp[s] += (1-0.5*S_dig[s])*GuoF(i,j,k,s,v[i,j,k])
for s in range(19):
f[i,j,k,s] = 0
for l in range(19):
#f=-M^-1*S(m-meq)
f[i,j,k,s] += inv_M[s,l]*m_temp[l]
periodic_index(i) set the bounary index with periodic bounary condition
@ti.func
def periodic_index(i):
#inner index
iout = i
#x-left
if i[0]<0: iout[0] = nx-1
#x-right
if i[0]>nx-1: iout[0] = 0
#y-left
if i[1]<0: iout[1] = ny-1
#y-right
if i[1]>ny-1: iout[1] = 0
#z-left
if i[2]<0: iout[2] = nz-1
#z-right
if i[2]>nz-1: iout[2] = 0
return iout
streaming1() defines the streaming process of denisty distibution function
@ti.kernel
def streaming1():
for i in ti.grouped(rho):
#if it is fluid
if (solid[i] == 0):
for s in range(19):
#the neighbour index
ip = periodic_index(i+e[s])
#if neighbour index is fluid just streaming
if (solid[ip]==0):
F[ip,s] = f[i,s]
#if neighbour index is solid just bounce back
else:
F[i,LR[s]] = f[i,s]
#print(i, ip, "@@@")
streaming2() a simple streaming process without consideration of solid and boundary
@ti.kernel
def streaming2():
for i in ti.grouped(rho):
for s in range(19):
f[i,s] = F[i,s]
Boudary_condition() define the bounary condition of fixed pressure and fixed velocity
@ti.kernel
def Boundary_condition():
#pressure-boundary condtion x-left
if ti.static(bc_x_left==1):
for j,k in ti.ndrange((0,ny),(0,nz)):
if (solid[0,j,k]==0):
for s in range(19):
#if boundary is fluid but the neighbour is solid
#equilibrium density distribution function is calculated based on the neighbour velocity
if (solid[1,j,k]>0):
F[0,j,k,s]=feq(s, rho_bcxl, v[1,j,k])
#if boundary is fluid and the neighbour is also fluid
#equilibrium density distribution function is calculated based on the boundary velocity
else:
F[0,j,k,s]=feq(s, rho_bcxl, v[0,j,k])
#velocity-boundary conditon x-left
if ti.static(bc_x_left==2):
for j,k in ti.ndrange((0,ny),(0,nz)):
if (solid[0,j,k]==0):
for s in range(19):
#calculate density distribution fucntion based on equilibrium part and non-equilibrium part
F[0,j,k,s]=feq(LR[s], 1.0, bc_vel_x_left[None])-F[0,j,k,LR[s]]+feq(s,1.0,bc_vel_x_left[None]) #!!!!!!change velocity in feq into vector
#pressure boundary condition x-right similar to x-left
if ti.static(bc_x_right==1):
for j,k in ti.ndrange((0,ny),(0,nz)):
if (solid[nx-1,j,k]==0):
for s in range(19):
if (solid[nx-2,j,k]>0):
F[nx-1,j,k,s]=feq(s, rho_bcxr, v[nx-2,j,k])
else:
F[nx-1,j,k,s]=feq(s, rho_bcxr, v[nx-1,j,k])
#velocity booundary condition x-right similar to x-left
if ti.static(bc_x_right==2):
for j,k in ti.ndrange((0,ny),(0,nz)):
if (solid[nx-1,j,k]==0):
for s in range(19):
F[nx-1,j,k,s]=feq(LR[s], 1.0, bc_vel_x_right[None])-F[nx-1,j,k,LR[s]]+feq(s,1.0,bc_vel_x_right[None]) #!!!!!!change velocity in feq into vector
streaming3() calculate the macroscopic variable
@ti.kernel
def streaming3():
for i in ti.grouped(rho):
#if it is fluid calculate density and velocity based on density distribution function
if (solid[i]==0):
rho[i] = 0
v[i] = ti.Vector([0,0,0])
for s in range(19):
f[i,s] = F[i,s]
rho[i] += f[i,s]
v[i] += e_f[s]*f[i,s]
v[i] /= rho[i]
v[i] += (ext_f[None]/2)/rho[i]
# if it is solid set denisty equals one and velocity equals zero
else:
rho[i] = 1.0
v[i] = ti.Vector([0,0,0])
At the end of the file do the actual simulation and export the data
#define some time varible
time_init = time.time()
time_now = time.time()
time_pre = time.time()
dt_count = 0
#import the solid flag data
#solid_np = init_geo('./BC.dat')
solid_np = init_geo('./img_ftb131.txt')
solid.from_numpy(solid_np)
# do the initialization
static_init()
init()
# do the actual simulation
for iter in range(50000+1):
colission()
streaming1()
Boundary_condition()
#streaming2()
streaming3()
# calculate every 1000 time step
if (iter%1000==0):
time_pre = time_now
time_now = time.time()
#calculate the time difference between now and previous time step
diff_time = int(time_now-time_pre)
#calculate the time difference between now and the initial time
elap_time = int(time_now-time_init)
#divmod function return the quotient and the remainder
#so that h_diff,m_diff and s_diff represent the hour, minute and second. the same as the h_elap,m_elap and s_elap
m_diff, s_diff = divmod(diff_time, 60)
h_diff, m_diff = divmod(m_diff, 60)
m_elap, s_elap = divmod(elap_time, 60)
h_elap, m_elap = divmod(m_elap, 60)
print('----------Time between two outputs is %dh %dm %ds; elapsed time is %dh %dm %ds----------------------' %(h_diff, m_diff, s_diff,h_elap,m_elap,s_elap))
print('The %dth iteration, Max Force = %f, force_scale = %f\n\n ' %(iter, 10.0, 10.0))
#export every 1000 timestep to vtk with x,y,z coordinate and solid,density and velocity variable
if (iter%10000==0):
gridToVTK(
"./structured"+str(iter),
x,
y,
z,
#cellData={"pressure": pressure},
pointData={ "Solid": np.ascontiguousarray(solid.to_numpy()),
"rho": np.ascontiguousarray(rho.to_numpy()),
"velocity": (np.ascontiguousarray(v.to_numpy()[:,:,:,0]), np.ascontiguousarray(v.to_numpy()[:,:,:,1]),np.ascontiguousarray(v.to_numpy()[:,:,:,2]))
}
)
# ti.sync()
# ti.profiler.print_kernel_profiler_info()
#print the profiler information of every kernel and task of taichi in this file
ti.profiler.print_scoped_profiler_info()