python进行OpenFOAM流场后处理

落花风

@田畔的风膜大佬！研究了半晚上。几点小疑问，向您请教：

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这条线。

田畔的风

@尚善若水画单纯的等值线可以参考 https://matplotlib.org/3.1.1/gallery/images_contours_and_fields/tricontour_smooth_delaunay.html#sphx-glr-gallery-images-contours-and-fields-tricontour-smooth-delaunay-py

尚善若水

@田畔的风好的，谢谢。

李东岳

fangyuanaza

@田畔的风这个流场图真的很漂亮。请教一下老师，我的网格是O型+H型，Fortran代码算的，能导入paraView。想请教一下，这种情况如何导入到python画图呢？如果直接用contourf 或者tricontourf都会导致网格变形。比如在paraview中为：
Screen Shot 2023-01-08 at 17.25.37.png
在python中直接
ax1.tricontourf(x,y, u,linewidth=0.25, levels=levels,cmap='coolwarm',vmin=0, vmax=1),显示出来为：
Screen Shot 2023-01-08 at 17.32.09.png
造成的原因可能是因为您说的Matplotlib只能处理三角形的网格。想请问在不清楚网格具体信息的情况下，能不能从paraview导出成三角网格？尝试过您代码中的Slice，但由于本身是2D，不能切片。请问有什么建议么？

田畔的风

我不太清楚你表述的网格扭曲是什么意思，如果只是纵横比不正确的话，可以使用set_aspect进行调整。比如我上面代码段中的ax.set_aspect('equal')，表示x轴和y轴始终保持1:1的比例。

fangyuanaza

@田畔的风感谢您的回复～

不是纵横比的问题，调整纵横比之后：
Screen Shot 2023-01-08 at 18.57.30.png

而是因为计算域边界是曲线，当调用tricontourf时，在一个矩形域内进行了插值，导致边界变形，成了直线。同理，用contourf时，首先需要对mesh进行插值，也会有同样的问题。

这个问题已经解决了，对于计算域的两个block, 分开处理。先reshapre到对应的block网格，然后用contourf来画。可以得到物理量的云图。

yuhxFoamer

@田畔的风
感谢老师分享非常实用的代码，本人也是初次使用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

田畔的风

@yuhxfoamer

首先pathlib这个模块在最终的代码中并没有被使用到，所以你可以将这个模块和对应的代码移除，这并不会影响代码运行结果。
然后谈谈你这个报错。readData.py这个文件是用pvpython执行的，这个python解释器是paraview自带的，对应的第三方库也是在安装paraview时配置好的，所以错误来源于你的paraview附带的pvpython解释器（可能是安装问题，可能是系统环境问题，也可能是软件版本问题，请自行排查），和系统python解释器无关。你使用pip和conda只是在安装你系统python解释器的第三方库。
尽管最终版本这个模块被弃用，不过这里再补充一下：pathlib在本代码中的用途是创建空文件夹（pathlib.Path），在上面的讨论中也提到了，因此可以使用os.mkdir替代。

知行

@落花风老师您好,请问下您是怎么安装的 paraview包呢,我在pypi里面没有发现这个包唉.非常感谢您.

CFD中文网

python进行OpenFOAM流场后处理