CFD中文网

    CFD中文网

    • 登录
    • 搜索
    • 最新

    python进行OpenFOAM流场后处理

    OpenFOAM
    7
    30
    1670
    正在加载更多帖子
    • 从旧到新
    • 从新到旧
    • 最多赞同
    回复
    • 在新帖中回复
    登录后回复
    此主题已被删除。只有拥有主题管理权限的用户可以查看。
    • F
      fangyuanaza 最后由 编辑

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

      李东岳 1 条回复 最后回复 回复 引用
      • 李东岳
        李东岳 管理员 @fangyuanaza 最后由 编辑

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

        CFD高性能服务器 http://dyfluid.com/servers.html

        F 1 条回复 最后回复 回复 引用
        • C
          chszkc 最后由 编辑

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

          F 落 2 条回复 最后回复 回复 引用
          • F
            fangyuanaza @chszkc 最后由 编辑

            @chszkc 感谢提供这个资源,试一试~

            1 条回复 最后回复 回复 引用
            • F
              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的。

              1 条回复 最后回复 回复 引用
              • 落
                落花风 @chszkc 最后由 编辑

                @chszkc 感谢分享,在尝试使用这个包的过程中遇到点小麻烦,向您请教。
                首先,我的计算域是这样的:19749d8e-3653-48bb-b6b4-a8946a742807.png

                我尝试画出温度等值线图,代码如下

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

                结果却是b3657482-2913-4d9f-aa47-59bbe9ee69ba.png

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

                C 1 条回复 最后回复 回复 引用
                • C
                  chszkc @落花风 最后由 chszkc 编辑

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

                  1 条回复 最后回复 回复 引用
                  • 田畔的风
                    田畔的风 讲师 最后由 编辑

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

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

                    QQ20220825-224603@2x.png
                    QQ20220825-224507@2x.png
                    QQ20220825-224633@2x.png

                    李东岳 Y 2 条回复 最后回复 回复 引用
                    • 田畔的风
                      田畔的风 讲师 最后由 编辑

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

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

                      落 1 条回复 最后回复 回复 引用
                      • 李东岳
                        李东岳 管理员 @田畔的风 最后由 编辑

                        @田畔的风 大佬后面3个图是什么出的,很好看,不像是apraview的风格

                        CFD高性能服务器 http://dyfluid.com/servers.html

                        田畔的风 1 条回复 最后回复 回复 引用
                        • 田畔的风
                          田畔的风 讲师 @李东岳 最后由 编辑

                          @李东岳 就是下面两段python代码绘制的,用的是matplotlib库。

                          1 条回复 最后回复 回复 引用
                          • 落
                            落花风 @田畔的风 最后由 编辑

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

                            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 条回复 最后回复 回复 引用
                            • 田畔的风
                              田畔的风 讲师 @落花风 最后由 田畔的风 编辑

                              @落花风

                              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自己的三角网格。

                              1. 如果你不关心切片后的几何拓扑,可以在slice函数中设置Triangulatetheslice=True。这个功能和paraview图形操作里是一样的,这样出来的切片就直接是三角网格了。就不需要进行上述的额外的三角化操作了。
                              落 F 2 条回复 最后回复 回复 引用
                              • 落
                                落花风 @田畔的风 最后由 编辑

                                @田畔的风 感谢,是我看的不仔细了。还有两点疑问:

                                     data = paraview.servermanager.Fetch(casefoam)
                                     data = dsa.WrapDataObject(data)
                                

                                这个data变量并没有 被用到,注释掉也没有影响(我是这样),data起什么作用呢?

                                Path("slice").mkdir(parents=True, exist_ok=True)
                                

                                建这个空文件夹是?

                                落 田畔的风 2 条回复 最后回复 回复 引用
                                • 落
                                  落花风 @落花风 最后由 编辑

                                  @落花风 放两张根据@田畔的风 大佬程序出的图,真的超棒。
                                  test3.png
                                  test4.svg

                                  1 条回复 最后回复 回复 引用
                                  • 李东岳
                                    李东岳 管理员 最后由 编辑

                                    太厉害了!一个教的细致一个学的细致。一个大作就出来了。

                                    CFD高性能服务器 http://dyfluid.com/servers.html

                                    1 条回复 最后回复 回复 引用
                                    • 田畔的风
                                      田畔的风 讲师 @落花风 最后由 编辑

                                      @落花风 代码只是一个早期雏形,很多调试功能最后并没有用到,我也懒得删干净了。不必理会。

                                      尚 1 条回复 最后回复 回复 引用
                                      • 尚
                                        尚善若水 @田畔的风 最后由 编辑

                                        @田畔的风 大佬,我不太懂python,正在照猫画虎修改您的代码。我想知道,b=10是什么意思呀?有没有办法默认读取latestTime结果。还有就是,如果我想在云图上加上另一个量的等值线,请问可以讲下如何操作吗?

                                        田畔的风 1 条回复 最后回复 回复 引用
                                        • 田畔的风
                                          田畔的风 讲师 @尚善若水 最后由 田畔的风 编辑

                                          @尚善若水

                                          其实是这么写的,就是指定读取的时间,后面改着改着就忘了

                                          b = 10
                                          
                                          os.system('pvpython readData.py -time_ %f' % b)
                                          

                                          你想读取最后一步也很简单。注意第一段代码中的这一行:

                                          times = reader.TimestepValues
                                          

                                          这一行代码返回了一个包含所有可读时间的List,List的最后一个元素就是latestTime。因此接下来使用这个元素更新pipeline即可:

                                          UpdatePipeline(time=times[-1])
                                          
                                          尚 1 条回复 最后回复 回复 引用
                                          • 尚
                                            尚善若水 @田畔的风 最后由 编辑

                                            @田畔的风 明白了,我看着应该是指定了特定时间,但是我不太确定。目前我网格数据4千多万,台式机paraview直接打不开了,正好看到了您的回答。请问我的第二个问题,添加等值线的可以指导下吗?谢谢。

                                            田畔的风 1 条回复 最后回复 回复 引用
                                            • 田畔的风
                                              田畔的风 讲师 @尚善若水 最后由 编辑

                                              @尚善若水 没能理解你的第二个问题

                                              尚 1 条回复 最后回复 回复 引用
                                              • 尚
                                                尚善若水 @田畔的风 最后由 编辑

                                                @田畔的风 比如我画出的了OH云图,但是我同时想在其上面绘制当量混合分数Zmix=constant这条线。

                                                田畔的风 1 条回复 最后回复 回复 引用
                                                • 田畔的风
                                                  田畔的风 讲师 @尚善若水 最后由 编辑

                                                  @尚善若水 画单纯的等值线可以参考 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

                                                  尚 1 条回复 最后回复 回复 引用
                                                  • 尚
                                                    尚善若水 @田畔的风 最后由 编辑

                                                    @田畔的风 好的,谢谢。

                                                    1 条回复 最后回复 回复 引用
                                                    • 李东岳
                                                      李东岳 管理员 最后由 编辑

                                                      :146: :146: :146: :146:

                                                      CFD高性能服务器 http://dyfluid.com/servers.html

                                                      1 条回复 最后回复 回复 引用
                                                      • F
                                                        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,不能切片。请问有什么建议么?

                                                        1 条回复 最后回复 回复 引用
                                                        • 田畔的风
                                                          田畔的风 讲师 最后由 编辑

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

                                                          F 1 条回复 最后回复 回复 引用
                                                          • F
                                                            fangyuanaza @田畔的风 最后由 编辑

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

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

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

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

                                                            1 条回复 最后回复 回复 引用
                                                            • Y
                                                              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
                                                              
                                                              田畔的风 1 条回复 最后回复 回复 引用
                                                              • 田畔的风
                                                                田畔的风 讲师 @yuhxFoamer 最后由 田畔的风 编辑

                                                                @yuhxfoamer

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