使用Python实现LBM(格子法)方腔驱动流

906

收藏 2019-12-14

<div class="pinggu_markdown">
<div class="pinggu_markdown__html"><h4 id="使用python实现lbm格子法方腔驱动流">使用Python实现LBM(格子法)方腔驱动流</h4>
<hr>
<ol>
<li>
<h5 id="python的不足：">Python的不足：</h5>
Python的最大优势也可能是它最大的弱点：它的灵活性和无类型的高级语法可能导致数据和计算密集型程序的性能不佳。—— 动态类型化解释语言
</li>
<li>
<h5 id="什么是-numba-">什么是 numba :</h5>
Numba，一个来自Anaconda的Python编译器，可以编译Python代码，以便在支持CUDA的GPU或多核CPU上执行。由于Python通常不是编译语言，您可能想知道为什么要使用Python编译器。答案当然是运行本机编译代码比运行动态解释代码快许多倍。 Numba允许您为Python函数指定类型签名，它可以在运行时进行编译（这是 “即时”或JIT编译）。 Numba动态编译代码的能力意味着您不会放弃Python的灵活性。这是向高效率编程和高性能计算提供理想组合的重要一步。
使用Numba，现在可以编写标准的Python函数并在支持CUDA的GPU上运行它们。 Numba专为面向阵列的计算任务而设计，就像广泛使用的NumPy库一样。面向阵列的计算任务中的数据并行性非常适合GPU等加速器。 Numba了解NumPy数组类型，并使用它们生成有效的编译代码，以便在GPU或多核CPU上执行。所需的编程工作可以像添加函数装饰器一样简单，以指示Numba为GPU编译。例如，以下代码中的@vectorize装饰器在运行时生成标量函数Add的已编译矢量化版本，以便它可用于在GPU上并行处理数据数组 
An example:
<pre class=" language-python"><code class="prism language-python"># -*- coding: utf-8 -*-
"""
Created on Thu Dec 5 14:00:59 2019

@author: Franz
"""
import numpy as np
from numba import vectorize

@vectorize(['float32(float32, float32)'], target='cpu')
def Add(a, b):
return a + b

# Initialize arrays
N = 100000
A = np.ones(N, dtype=np.float32)
B = np.ones(A.shape, dtype=A.dtype)
C = np.empty_like(A, dtype=A.dtype)

# Add arrays on GPU
C = Add(A, B)
</code></pre>
要在CPU上编译和运行相同的函数，我们只需将目标更改为“cpu”，从而在CPU上编译的矢量化C代码级别上产生性能,灵活性大大提高。
</li>
<li>
numpy的<code>loadtxt()</code>函数的用法 
<code>numpy.loadtxt(fname, dtype=, comments='#', delimiter=None, converters=None, skiprows=0, usecols=None, unpack=False, ndmin=0)</code>

<table>
<thead>
<tr>
<th align="center">参数</th>
<th align="center">作用</th>
</tr>
</thead>
<tbody>
<tr>
<td align="center">fname</td>
<td align="center">被读取的文件名（相对地址与绝对地址<a href="#fn1" id="fnref1">1</a>）</td>
</tr>
<tr>
<td align="center">comments</td>
<td align="center">跳过以指定元素开头的行（不读取）</td>
</tr>
<tr>
<td align="center">delimiter</td>
<td align="center">指定读取文件中的数据分割符</td>
</tr>
<tr>
<td align="center">converters<a href="#fn2" id="fnref2">2</a></td>
<td align="center">对读取数据进行预处理</td>
</tr>
<tr>
<td align="center">skiprows</td>
<td align="center">选择跳过的行数</td>
</tr>
<tr>
<td align="center">usecols</td>
<td align="center">指定读取列</td>
</tr>
<tr>
<td align="center">unpack</td>
<td align="center">是否将数据进行向量化输出<a href="#fn3" id="fnref3">3</a></td>
</tr>
<tr>
<td align="center">encoding</td>
<td align="center">对数据进行预编码<a href="#fn4" id="fnref4">4</a></td>
</tr>
</tbody>
</table></li>
<li>
用matplotlib.pyplot的quiver()函数绘制矢量图 
调用方式
<pre class=" language-python"><code class="prism language-python">quiver(U, V, **kw)
quiver(U, V, C, **kw)
quiver(X, Y, U, V, **kw)
quiver(X, Y, U, V, C, **kw)
</code></pre>
U、V是箭头数据（data），X、Y是箭头的位置，C是箭头的颜色。这些参数可以是一维或二维的数组或序列。 
如果X、Y不存在（absent），他们将作为均匀网格被生成。如果U和V是2维的数组，X和Y是1维数组，并且len(X)和len(Y)与U的列（column）和行（row）纬度相匹配（match），那么X和Y将使用numpy.meshgrid()——用于产生一个矩阵，具体可参考：meshgrid使用方法——进行扩展。
默认设置会自动将箭头的长度缩放到合理的大小。要改变箭头的长度请参看 scale 和scale_units两个关键字参数（kwargs:关键字参数，参看文章最后有关键字参数与可变参数的区别）
默认值给出一个稍微后掠（swept-back）的箭头；若要使箭头头部呈三角状，则要确保headaxislength与headlength相同。若要使箭头更尖锐（more pointed），则应减小headwidth或者增大headlength和headaxislength。若要使箭头头部相对于箭杆（shaft）更小一些，则应将所有头部参数都减小（scale down）。你最好不要单独留下minshaft（You will probably do best to leave minshaft alone.） 
相关详细用法可以参照<a href="https://blog.csdn.net/liuchengzimozigreat/article/details/84566650">Blog</a>,以及官方的<a href="https://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.quiver">matplotlib参考文档</a>
</li>
<li>
通过上述简单的加速，并通过对顶盖驱动流改写成Python程序，通过对其使用Numpy和Numba进行加速，最后使用quiver（）函数绘制流场图：
<pre class=" language-python"><code class="prism language-python"># -*- coding: utf-8 -*-
"""
Created on Fri Dec 6 15:17:10 2019

@author: Franz
"""
import numpy as np
import numba as nb
import matplotlib.pyplot as plt

plt.rcParams['figure.figsize'] = (5.0, 5.0)
plt.rcParams['figure.dpi'] = 100 #分辨率
# function
# 1st:equilibrum function
@nb.jit()
def Feq(k,rho,u):
eu = e[k,0]*u[0] + e[k,1]*u[1];
uv = u[0]**2 + u[1]**2;
feq = rho*w[k]*(1.0+3.0*eu+4.5*eu*eu-1.5*uv);
return feq

# 2nd:evolution: including move and collide,calculate micro-parameter
@nb.jit()
def Evolution(f,rho,u):
F = np.zeros((NX+1,NY+1,Q));
for i in range(1,NX):
 for j in range(1,NY):
 for k in range(Q):
 ip = i - e[k,0];
 jp = j - e[k,1];
 F[i,j,k] = f[ip,jp,k] + (Feq(k,rho[ip,jp],u[ip,jp])-f[ip,jp,k])/tau_f;

u0 = np.empty((NX+1,NY+1,2));
for i in range(1,NX):
 for j in range(1,NY):
 u0[i,j,0] = u[i,j,0];
 u0[i,j,1] = u[i,j,1];
 rho[i,j] = 0;
 u[i,j,0] = 0;
 u[i,j,1] = 0;
 for k in range(Q):
 f[i,j,k] = F[i,j,k];
 rho[i,j] += f[i,j,k];
 u[i,j,0] += e[k,0]*f[i,j,k];
 u[i,j,1] += e[k,1]*f[i,j,k];
 u[i,j,0] /= rho[i,j];
 u[i,j,1] /= rho[i,j];

# left & right marging
for j in range(1,NY):
 for k in range(Q):
 rho[NX,j] = rho[NX-1,j];
 f[NX,j,k]=Feq(k,rho[NX,j],u[NX,j])+f[NX-1,j,k]-Feq(k,rho[NX-1,j],u[NX-1,j]);
 rho[0,j]=rho[1,j];
 f[0,j,k]=Feq(k,rho[0,j],u[0,j])+f[1,j,k]-Feq(k,rho[1,j],u[1,j]);

# top & bottom margin
for i in range(NX+1):
 for k in range(Q):
 rho[i,0] = rho[i,1];
 f[i,0,k]=Feq(k,rho[i,0],u[i,0])+f[i,1,k]-Feq(k,rho[i,1],u[i,1]);
 rho[i,NY]=rho[i,NY-1];
 u[i,NY,0]=U;
 f[i,NY,k]=Feq(k,rho[i,NY],u[i,NY])+f[i,NY-1,k]-Feq(k,rho[i,NY-1],u[i,NY-1]);

return f,u,u0

@nb.jit()
def Error(u,u0):
temp1 = 0
temp2 = 0
for i in range(1,NX):
 for j in range(1,NY):
 temp1 += (u[i,j,0]-u0[i,j,0])**2+(u[i,j,1]-u0[i,j,1])**2;
 temp2 += u[i,j,0]**2+u[i,j,1]**2;
 error = np.sqrt(temp1)/np.sqrt(temp2+1e-30);
return error


global Q,NX,NY,U;
Q = 9;
NX = 256;
NY = 256;
U = 0.1;

global e,w,tau_f;
e = np.array([[0,0],[1,0],[0,1],[-1,0],[0,-1],[1,1],[-1,1],[-1,-1],[1,-1]],dtype=int);
w = np.array([4.0/9,1.0/9,1.0/9,1.0/9,1.0/9,1.0/36,1.0/36,1.0/36,1.0/36]);

# main
# init
dx = 1.0;
dy = 1.0;
Lx = dx * NX;
Ly = dy * NY;
dt = dx;
c = dx/dt;
rho0 = 1.0;
Re = 1000;
niu = U * Lx / Re;
tau_f = 3.0*niu + 0.5;

u = np.zeros((NX+1,NY+1,2),dtype='double');
rho = np.ones((NX+1,NY+1),dtype='double')*rho0;
f = np.empty((NX+1,NY+1,Q),dtype='double');
for i in range(NX+1):
u[i,NY,0]=U;
for j in range(NY+1):
 for k in range(Q):
 f[i,j,k] = Feq(k,rho[i,j],u[i,j]);

n = 1; # calculate the times
while True:
f,u,u0 = Evolution(f,rho,u);
if n % 100 == 0:
 if Error(u,u0) < 1.0e-6:
 break
 if n >= 1000 & n % 1000 == 0:
 x,y = np.mgrid[0:257:257J,0:257:257J];
 ux = u[:,:,0];
 uy = u[:,:,1];
 M = np.hypot(ux, uy)
 a = plt.quiver(ux[::9,::9].T,uy[::9,::9].T,M[::9,::9])
 plt.show()
 print('The run times is %d'%(n));
n += 1;
</code></pre>
</li>
</ol>
<hr class="footnotes-sep">
<section class="footnotes">
<ol class="footnotes-list">
<li id="fn1" class="footnote-item">绝对路径和相对路径的区别及其转换可以查看<a href="https://www.cnblogs.com/wangyanyan/p/7440685.html">Blog</a> <a href="#fnref1" class="footnote-backref">↩︎</a>
</li>
<li id="fn2" class="footnote-item">converts对格式的控制类似<code>converters={1:_is_num}</code>,其中1代表控制输出第一列,，详细参考<a href="https://www.cnblogs.com/maoguy/p/6894591.html">Blog</a> <a href="#fnref2" class="footnote-backref">↩︎</a>
</li>
<li id="fn3" class="footnote-item"><code>unpack</code>：选择是否将数据向量输出，默认是False，即将数据逐行输出，当设置为True时，数据将逐列输出 <a href="#fnref3" class="footnote-backref">↩︎</a>
</li>
<li id="fn4" class="footnote-item"><code>encoding</code>决定读取文件时使用的编码方式,可以参照<a href="https://blog.csdn.net/lwhello/article/details/80271633">Blog</a> <a href="#fnref4" class="footnote-backref">↩︎</a>
</li>
</ol>
</section>
</div>
</div>

扫码加我拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

栏目导航

扫码加我 拉你入群

分享

扫码加好友，拉您进群

扫码加我拉你入群