python进行OpenFOAM流场后处理
-
请教各位老师,在做几组流场图对比的时候,不用软件导出来的图,而是导出数据,用python进行绘制,改如何做呢?比如tecplot形式生成的数据。
-
openfoam的结果都存储在时间步文件夹。你所需要的大体上就是对这些文件序列做处理了吧?但是如果你要抽取某条线上的值,大概率还需要paraview
-
https://fluidfoam.readthedocs.io/en/latest/
我一直都是用这个包来读OpenFOAM cell的数据来做后处理的,非常好用 -
@chszkc 感谢提供这个资源,试一试~
-
@李东岳 是的,我当时不清楚该怎么做来用python处理整个流场数据,最后是用paraview导入tecplot数据,然后导出成csv等方便读取数据的格式,用python进行处理。因为tecplot数据结构是分块的,而csv等格式是坐标点对应物理量的形式,方便读取。
第二个可能对大家有用的地方在于:有时候tecplot数据导入paraview会一直卡顿,无法完全显示,需要删掉一些符号:
• X, Y and Z array headers MUST NOT have other characters like [mm]. It must be simply “x”, “y” and “z”.Paraview also accepts this: “Z, mm” and this “x/c” • This parameter is not allowed by Paraview (it must be deleted from the file): “ZONETYPE=…” • “STRANDID” and “SOLUTION TIME” are currently unsupported by Paraview but they can be left in the file (it gives only a warning message) • Comments are allowed by using “#” (like Python)我当时是删除了ZONETYPE, 才能导入paraview的。
-
@chszkc 感谢分享,在尝试使用这个包的过程中遇到点小麻烦,向您请教。
首先,我的计算域是这样的:我尝试画出温度等值线图,代码如下
#%% from fluidfoam import readmesh sol = './T293/massFlowRate0.412' x, y, z = readmesh(sol) from fluidfoam import readscalar timename = 'latestTime' T = readscalar(sol, timename, 'T') # %% import numpy as np import matplotlib.pyplot as plt # Define plot parameters fig, ax = plt.subplots(figsize=(7, 3), dpi=100) plt.rcParams.update({'font.size': 10}) plt.xlabel('x') plt.ylabel('y') # Plots the contour of sediment concentration levels = np.arange(30, 300, 10) plt.tricontourf(x, y, T, cmap=plt.cm.Reds, levels=levels)
结果却是
请问您有类似经历吗
先谢谢您 -
@落花风 抱歉我暂时没用python画过这么复杂几何+非结构网格的云图....我的算例是基于均匀网格的所以相对比较好处理
-
比较方便的方法是联合paraview的pvpython和你自己的python进行处理:
- 使用pvpython做切片等一系列操作,导出一个临时的数据文件。
- 然后使用你自己的python,利用meshio包读取之前的数据文件,然后用matplotlib绘制云图(仅支持三角形,多边形可以自己对单元做个三角化处理下)
-
代码分两个文件,一个读数据,一个绘制。
readData.py
from paraview.simple import * from pathlib import Path import sys from paraview.vtk.numpy_interface import dataset_adapter as dsa # CaseType: 0 = decomposed case, 1 = reconstructed case. reader = OpenFOAMReader(FileName="PEE.foam", CaseType=0) if reader: print("Successful read.") else: print("Failed read.") time = reader.TimestepValues # print(time) argList = sys.argv # print(argList) try: timeIndex = argList.index('-time_') + 1 timeI = float(argList[timeIndex]) except ValueError: print('ValueError: undefined time parameter.') exit() UpdatePipeline(time=timeI) data = paraview.servermanager.Fetch(reader) data = dsa.WrapDataObject(data) slicer = Slice(Input=reader, SliceType="Plane", Triangulatetheslice=False) slicer.SliceType.Origin = [0, 5.5, 0] slicer.SliceType.Normal = [0, 1, 0] Path("slice").mkdir(parents=True, exist_ok=True) SaveData("/Volumes/RAM Disk/slicePlane_%.2f.e" % timeI, slicer) print('Successful write.')
plotData.py
import meshio import numpy as np import matplotlib.pyplot as plt import matplotlib.tri as tri import os #from matplotlib import ticker, cm plt.rcParams.update({ # "text.usetex": True, "font.family": "serif", "font.serif": ['SimSun'], # "axes.unicode_minus": False, "mathtext.fontset": "stix", "font.size": 16, "figure.figsize": (12, 5)}) ticklabel_style = {"fontname": "Times New Roman", "fontsize": 14} b = 10 os.system('pvpython readData.py -time_ 10') mesh = meshio.read('/Volumes/RAM Disk/slicePlane_%.2f.e' % b) os.system('rm -r "/Volumes/RAM Disk/slicePlane_%.2f.e"' % b) logData = np.where(mesh.point_data['driftDensity'] < 1e-12, -12., np.log10(mesh.point_data['driftDensity'])) fig, ax = plt.subplots() ax.set_aspect('equal') cell_0 = mesh.cells[0].data[:, :3] print(mesh.cells[0].data[:, 2:]) print(mesh.cells[0].data[:, 0]) cell_1 = np.concatenate( (mesh.cells[0].data[:, 2:], mesh.cells[0].data[:, [0]]), axis=1) cell = np.concatenate((cell_0, cell_1), axis=0) triang = tri.Triangulation(mesh.points[:, 0], mesh.points[:, 2], cell) #triang = tri.Triangulation(mesh.points[:,0], mesh.points[:,2]) tpc = ax.tricontourf(triang, logData, cmap='jet', levels=16) #ax.triplot(triang, 'k-', linewidth=0.5) for cellI in cell: ax.plot(mesh.points[:, 0][cellI], mesh.points[:, 2][cellI],'k-', linewidth=0.5) cbar = fig.colorbar(tpc, format='$10^{%.1f}$', fraction=0.0172, pad=0.06) cbar.ax.tick_params(labelsize=14) plt.savefig('test.jpg', dpi=400) plt.show() # print(mesh.cell_data)
注意我这里的存储路径
/Volumes/RAM Disk/
是macOS风格的,Linux下面可以把临时文件存在/dev/shm/
,这样可以回避外部存储,直接在内存中实现IO。 -
@田畔的风 大佬后面3个图是什么出的,很好看,不像是apraview的风格
-
@李东岳 就是下面两段python代码绘制的,用的是matplotlib库。
-
@田畔的风 膜大佬!研究了半晚上。几点小疑问,向您请教:
reader = OpenFOAMReader(FileName="PEE.foam", CaseType=0)
- caseType=0 是控制什么的?这句导致我的case读不到任何一个时间步。删除之后可以读到。
logData = np.where( mesh.point_data['driftDensity'] < 1e-12, -12.0, np.log10(mesh.point_data['driftDensity']), )
- 这句是控制什么的呢?我print(mesh.dict),发现没有‘driftdensity’
- 我用的是六面体网格,怎么划分成三角形,你的代码里有相关的吗?(我刚理解到第二个问题,后面还没看,肝不动了)
多谢分享,多谢多谢
-
-
caseType
参数的作用在上一行的注释中已经给出。caseType=0
指读写decompose
算例,即包含processor*
文件夹的并行算例。caseType=1
指读写reconstruct
算例,即串行或使用reconstructPar
重建后的算例。 -
driftDensity
是我自己写的某个求解器里的标量场的名称。在标准求解器里,你想读压力就是p
,想读速度就是U
,想读温度就是T
... -
我的这个算例用的也是六面体网格,通过一次
slice
操作转换成了四边形切片,然后对每个四边形单元进行三角化。请注意理解我的这段代码:
cell_0 = mesh.cells[0].data[:, :3] cell_1 = np.concatenate( (mesh.cells[0].data[:, 2:], mesh.cells[0].data[:, [0]]), axis=1) cell = np.concatenate((cell_0, cell_1), axis=0) triang = tri.Triangulation(mesh.points[:, 0], mesh.points[:, 2], cell)
就是取每个四边形的
0 1 2
和0 2 3
两组顶点,组成两个互不重叠的三角形,然后生成matplotlib自己的三角网格。- 如果你不关心切片后的几何拓扑,可以在slice函数中设置
Triangulatetheslice=True
。这个功能和paraview图形操作里是一样的,这样出来的切片就直接是三角网格了。就不需要进行上述的额外的三角化操作了。
-
-
@田畔的风 感谢,是我看的不仔细了。还有两点疑问:
data = paraview.servermanager.Fetch(casefoam) data = dsa.WrapDataObject(data)
这个data变量并没有 被用到,注释掉也没有影响(我是这样),data起什么作用呢?
Path("slice").mkdir(parents=True, exist_ok=True)
建这个空文件夹是?
-
@落花风 放两张根据@田畔的风 大佬程序出的图,真的超棒。
-
太厉害了!一个教的细致一个学的细致。一个大作就出来了。
-
@落花风 代码只是一个早期雏形,很多调试功能最后并没有用到,我也懒得删干净了。不必理会。
-
@田畔的风 大佬,我不太懂python,正在照猫画虎修改您的代码。我想知道,b=10是什么意思呀?有没有办法默认读取latestTime结果。还有就是,如果我想在云图上加上另一个量的等值线,请问可以讲下如何操作吗?
-
其实是这么写的,就是指定读取的时间,后面改着改着就忘了
b = 10 os.system('pvpython readData.py -time_ %f' % b)
你想读取最后一步也很简单。注意第一段代码中的这一行:
times = reader.TimestepValues
这一行代码返回了一个包含所有可读时间的List,List的最后一个元素就是latestTime。因此接下来使用这个元素更新pipeline即可:
UpdatePipeline(time=times[-1])
-
@田畔的风 明白了,我看着应该是指定了特定时间,但是我不太确定。目前我网格数据4千多万,台式机paraview直接打不开了,正好看到了您的回答。请问我的第二个问题,添加等值线的可以指导下吗?谢谢。
-
@尚善若水 没能理解你的第二个问题
-
@田畔的风 比如我画出的了OH云图,但是我同时想在其上面绘制当量混合分数Zmix=constant这条线。
-
-
@田畔的风 好的,谢谢。
-
-
@田畔的风 这个流场图真的很漂亮。请教一下老师,我的网格是O型+H型,Fortran代码算的,能导入paraView。想请教一下,这种情况如何导入到python画图呢?如果直接用contourf 或者tricontourf都会导致网格变形。比如在paraview中为:
在python中直接
ax1.tricontourf(x,y, u,linewidth=0.25, levels=levels,cmap='coolwarm',vmin=0, vmax=1),显示出来为:
造成的原因可能是因为您说的Matplotlib只能处理三角形的网格。想请问在不清楚网格具体信息的情况下,能不能从paraview导出成三角网格?尝试过您代码中的Slice,但由于本身是2D,不能切片。请问有什么建议么? -
我不太清楚你表述的网格扭曲是什么意思,如果只是纵横比不正确的话,可以使用
set_aspect
进行调整。比如我上面代码段中的ax.set_aspect('equal')
,表示x轴和y轴始终保持1:1的比例。 -
@田畔的风 感谢您的回复~
不是纵横比的问题,调整纵横比之后:
而是因为计算域边界是曲线,当调用tricontourf时,在一个矩形域内进行了插值,导致边界变形,成了直线。同理,用contourf时,首先需要对mesh进行插值,也会有同样的问题。
这个问题已经解决了,对于计算域的两个block, 分开处理。先reshapre到对应的block网格,然后用contourf来画。可以得到物理量的云图。
-
@田畔的风
感谢老师分享非常实用的代码,本人也是初次使用pvpython处理数据,请问老师一个简单的问题,当我运行您的程序时,出现了以下的错误,我想这个错误应该是非常常见的,但是目前我尝试的方法有:1 重新install pathlib(用pipe 和 conda都试过) 2. 更新python的版本,这两个方法都没有将此错误排除,求问老师有没有其他方法?Traceback (most recent call last): File "./readData.py", line 3, in <module> from pathlib import Path ImportError: No module named pathlib
-
- 首先pathlib这个模块在最终的代码中并没有被使用到,所以你可以将这个模块和对应的代码移除,这并不会影响代码运行结果。
- 然后谈谈你这个报错。
readData.py
这个文件是用pvpython
执行的,这个python解释器是paraview自带的,对应的第三方库也是在安装paraview时配置好的,所以错误来源于你的paraview附带的pvpython解释器(可能是安装问题,可能是系统环境问题,也可能是软件版本问题,请自行排查),和系统python解释器无关。你使用pip
和conda
只是在安装你系统python解释器的第三方库。 - 尽管最终版本这个模块被弃用,不过这里再补充一下:pathlib在本代码中的用途是创建空文件夹(
pathlib.Path
),在上面的讨论中也提到了,因此可以使用os.mkdir
替代。