python进行OpenFOAM流场后处理

[fangyuanaza]

请教各位老师，在做几组流场图对比的时候，不用软件导出来的图，而是导出数据，用python进行绘制，改如何做呢？比如tecplot形式生成的数据。

[李东岳]

openfoam的结果都存储在时间步文件夹。你所需要的大体上就是对这些文件序列做处理了吧？但是如果你要抽取某条线上的值，大概率还需要paraview

[chszkc]

https://fluidfoam.readthedocs.io/en/latest/
我一直都是用这个包来读OpenFOAM cell的数据来做后处理的，非常好用

[fangyuanaza]

@chszkc 感谢提供这个资源，试一试～

[fangyuanaza]

@李东岳 是的，我当时不清楚该怎么做来用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)

结果却是

请问您有类似经历吗
先谢谢您

[chszkc]

@落花风 抱歉我暂时没用python画过这么复杂几何+非结构网格的云图....我的算例是基于均匀网格的所以相对比较好处理

[田畔的风]

比较方便的方法是联合paraview的pvpython和你自己的python进行处理：

1. 使用pvpython做切片等一系列操作，导出一个临时的数据文件。
2. 然后使用你自己的python，利用meshio包读取之前的数据文件，然后用matplotlib绘制云图（仅支持三角形，多边形可以自己对单元做个三角化处理下）

[田畔的风]

代码分两个文件，一个读数据，一个绘制。

1. 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.')

2. 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’
• 我用的是六面体网格，怎么划分成三角形，你的代码里有相关的吗？（我刚理解到第二个问题，后面还没看，肝不动了）

多谢分享，多谢多谢

[田畔的风]

@落花风

1. caseType参数的作用在上一行的注释中已经给出。caseType=0指读写decompose算例，即包含processor*文件夹的并行算例。caseType=1指读写reconstruct算例，即串行或使用reconstructPar重建后的算例。

2. driftDensity是我自己写的某个求解器里的标量场的名称。在标准求解器里，你想读压力就是p，想读速度就是U，想读温度就是T...

3. 我的这个算例用的也是六面体网格，通过一次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自己的三角网格。

4. 如果你不关心切片后的几何拓扑，可以在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直接打不开了，正好看到了您的回答。请问我的第二个问题，添加等值线的可以指导下吗？谢谢。