Numba学习日记 —— 2019-12-5

Numba学习日记 —— 2019-12-5

1. Python的不足：

   Python的最大优势也可能是它最大的弱点：它的灵活性和无类型的高级语法可能导致数据和计算密集型程序的性能不佳。—— 动态类型化解释语言

2. 什么是 numba :

   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:

   # -*- 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)

   要在CPU上编译和运行相同的函数，我们只需将目标更改为“cpu”，从而在CPU上编译的矢量化C代码级别上产生性能,灵活性大大提高。

3. numpy的loadtxt()函数的用法
   numpy.loadtxt(fname, dtype=, comments=‘#‘, delimiter=None, converters=None, skiprows=0, usecols=None, unpack=False, ndmin=0)

   参数	作用
   fname	被读取的文件名（相对地址与绝对地址1）
   comments	跳过以指定元素开头的行（不读取）
   delimiter	指定读取文件中的数据分割符
   converters2	对读取数据进行预处理
   skiprows	选择跳过的行数
   usecols	指定读取列
   unpack	是否将数据进行向量化输出3
   encoding	对数据进行预编码4

4. 用matplotlib.pyplot的quiver()函数绘制矢量图
   调用方式
   python quiver(U, V, **kw) quiver(U, V, C, **kw) quiver(X, Y, U, V, **kw) quiver(X, Y, U, V, C, **kw)
   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.）
   相关详细用法可以参照Blog,以及官方的matplotlib参考文档

5. 通过上述简单的加速，并通过对顶盖驱动流改写成Python程序，通过对其使用Numpy和Numba进行加速，最后使用quiver（）函数绘制流场图：

   # -*- 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;

原文地址：https://www.cnblogs.com/lovely-bones/p/12001844.html