python输出方格_使用Python实现LBM(格子法)方腔驱动流
使用Python实现LBM(格子法)方腔驱动流
Python的不足:
Python的最大优势也可能是它最大的弱点:它的灵活性和无类型的高级语法可能导致数据和计算密集型程序的性能不佳。—— 动态类型化解释语言
什么是 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代码级别上产生性能,灵活性大大提高。
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是否将数据进行向量化输出3encoding对数据进行预编码4
用matplotlib.pyplot的quiver()函数绘制矢量图
调用方式
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参考文档
通过上述简单的加速,并通过对顶盖驱动流改写成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;
绝对路径和相对路径的区别及其转换可以查看Blog ↩︎
converts对格式的控制类似
converters={1:_is_num}
,其中1代表控制输出第一列,,详细参考Blog ↩︎unpack
:选择是否将数据向量输出,默认是False,即将数据逐行输出,当设置为True时,数据将逐列输出 ↩︎encoding
决定读取文件时使用的编码方式,可以参照Blog ↩︎
python输出方格_使用Python实现LBM(格子法)方腔驱动流相关推荐
- python输出运行时间_分析python程序运行时间的几种方法
最早见过手写的,类似于下面这种: 1 import datetime 2 3 def time_1(): 4 begin = datetime.datetime.now() 5 sum = 0 6 f ...
- 格子玻尔兹曼机(Lattice Boltzmann Method)系列4:LBM实例之方腔驱动流
问题描述 在一个正方形的空腔中,顶盖以U的速度被突然启动,导致产生了一个在边角处会有小的漩涡出现的流场. 左.右.下边界使用了回弹性边界条件,而上边界使用了上一篇文章中的Von Neumann速度边界 ...
- python 概率分布模型_使用python的概率模型进行公司估值
python 概率分布模型 Note from Towards Data Science's editors: While we allow independent authors to publis ...
- python 时间序列预测_使用Python进行动手时间序列预测
python 时间序列预测 Time series analysis is the endeavor of extracting meaningful summary and statistical ...
- python 打印类型_让Python输出更漂亮:PrettyPrinter
PrettyPrinter是Python 3.6 及以上版本中的一个功能强大.支持语法高亮.描述性的美化打印包.它使用了改进的Wadler-Leijen布局算法,和Haskell打印美化库中的pret ...
- python高斯求和_利用Python进行数据分析(3)- 列表、元组、字典、集合
本文主要是对Python的数据结构进行了一个总结,常见的数据结构包含:列表list.元组tuple.字典dict和集合set. image 索引 左边0开始,右边-1开始 通过index()函数查看索 ...
- python希腊字母字符串_#10 Python字符串
前言 通过上一节可知,Python6个序列的内置类型中,最常见的是列表和元组,但在Python中,最常用的数据类型却不是列表和元组,而是字符串.要想深入了解字符串,必须先掌握字符编码问题.因此本篇博文 ...
- python选择题题目_《Python程序设计》题库 - 选择题
一.基础知识 1 . Python 语言属于( ) . C A . 机器语言 B . 汇编语言 C .高级语言 D .科学计算语言 2 .下列选项中,不属于 Python 特点的是( ) . B A ...
- 用python实现点阵屏_用Python代码来绘制彭罗斯点阵的教程
这里是显示彭罗斯点阵的Python的脚本.是的,这是可以运行的有效Phython代码. 译注:彭罗斯点阵,物理学术语.上世纪70年代英国数学家彭罗斯第一次提出了这个概念,称为彭罗斯点阵(Pen-ros ...
最新文章
- 微软职位内部推荐-Software Development Engineer II
- object.create()
- sqlplus 如何连接到指定数据库,并创建用户与授权
- 人工智能助力资深内容营销人员
- x86服务器当虚拟化的存储,X86服务器虚拟化实施方案.doc
- JS中深浅拷贝 函数封装代码
- [react] 怎样将事件传递给子组件?
- 深入java核心_Java核心(五)深入理解BIO、NIO、AIO
- 用python实现时间的动态(动态时钟)+ 算出某年某月星期几的所有日期
- Qt笔记-添加Win10Pcap库获取网络适配器(MinGW编译器)
- 探究foreach对于迭代变量的封装性的研究
- 玩转SpringSession,重要知识点全面剖析!
- 【业界】百度NLP十年技术积累,最新发布5款产品,公布两大计划
- 遇见Flask-Script
- 计算机组成白中英答案,计算机组成原理白中英答案
- Java开发实战经典【Java基础】
- MAC协议微信自动生成带备注固定收款码免挂机监控协议回调
- Excel如何合并单元格
- EcmaScript 2022中的新特性
- 2021中国华录杯·算法大赛直通车!
热门文章
- Mali GPU OpenGL ES 应用性能优化--基本概念
- 【网络模拟】Network Emulator for Windows Toolkit
- QQ靓号申请器v1.2.0.0【源码】
- 大数据是如何进行分析的
- 骑马式瘦腿瑜伽 消除肌肉型小腿
- 制作IE Cab包的过程
- python编程实现将文本音频数据还原为wav语音文件
- WIN10+MATLAB2018b+STK11.6+MATLAB_Connectors1.0.11安装记录
- VirtuaNES.v0.97源码探究6 内存相关
- 【51nod_3202】子集和判断