ARTICLE DETAIL

资讯详情

深耕网站视觉设计与运营推广的一线实战洞察。

OpenFOAM二次开发教程(17):Python 前后处理(二)——pyvista/pandas 结果聚合与可视化

OpenFOAM二次开发教程(17):Python 前后处理(二)——pyvista/pandas 结果聚合与可视化 OpenFOAM二次开发教程17Python 前后处理二——pyvista/pandas 结果聚合与可视化版本与事实声明函数对象的产物组织与执行方式见官方文档站Function objects页postProcess工具无附加选项时执行 controlDict 中为所有时间目录列出的函数对象的说明见官方文档站postProcess页。postProcess是官方后处理工具官方 API 文档另有foamPostProcess应用二者用途相关、具体用法以官方文档为准。foamToVTK是 OpenFOAM 常用的结果转换工具用于导出为 VTK 格式供外部可视化工具读取具体命令行选项以官方文档与你本机-help为准。pyvista、pandas、matplotlib均为 PyPI 上的第三方库非 OpenFOAM 官方发布物foamlib亦来自 PyPI。文本对官方与第三方的区分严格标明。文中数值与坐标为示例不代表任何标准规定。一句话结论OpenFOAM 的后处理数据有两条出口——postProcessing/下的函数对象时间序列结构化、适合pandas聚合与场数据适合pyvista三维可视化把这两条出口接进 Python就能实现多算例自动汇总成表 → 生成对比图 → 出报告的完整闭环而唯一不可妥协的纪律是每个数字都能回溯到某个算例的某个文件。〇、本篇要解决的认知问题Q1OpenFOAM 的后处理数据到底有哪几种形态各自适合什么工具Q2postProcessing/下的文件是什么结构为什么它比解析日志可靠Q3怎么用pandas把多算例、多时刻的监测数据聚合成一张对比长表Q4怎么把场数据读进pyvista做三维可视化Q5出图出表时哪些纪律决定了这张图能不能被信任一、机制解析1.1 后处理数据的三种形态OpenFOAM 运行结束后数据以三种形态存在形态位置结构最适合时间序列函数对象输出postProcessing/函数对象名/时间/下的文件常为 CSV 或数据文件结构化、按时间排列pandas聚合、趋势对比、出图场数据时间目录/如500/U、500/pOpenFOAM 自有文本/二进制场格式三维可视化、空间分析日志文本log.求解器非结构化文本快速排查、残差轨迹不适合作为长期数据源关键判断本篇最重要的一条能用功能对象产出结构化数据就不要去解析日志。官方的函数对象第 12 篇会把监测数据写成结构化文件而日志格式随版本演进、随求解器不同而变化——基于日志的解析脚本是最脆弱的技术债。迁移动作如果你现在的后处理是grep日志 awk取数请把它改造成用probes/surfaceFieldValue等函数对象产出 CSV pandas读 CSV。这一步的收益是长期的。1.2 postProcessing/ 的组织结构与 postProcess 工具函数对象的产物一般落在postProcessing/下按函数对象名/时间/文件名分层组织。具体目录层级与文件命名随函数对象类型与版本不同——因此解析脚本应当先探测目录结构、再按模式匹配文件而不是硬编码一个路径。官方文档对postProcess工具的说明是无附加选项时它执行controlDict中为所有时间目录列出的函数对象。这带来一个非常重要的能力同一个函数对象配置既能在线在求解过程中边算边输出也能离线在求解结束后对已有时间目录补算。实践价值算例跑完了才发现忘了加某个监测不用重跑——用postProcess离线执行即可。这是把后处理从必须提前规划变成可以事后补的关键机制。1.3 pandas 聚合从一个算例一条曲线到一张对比表聚合的规范做法是长表long formatcase param_value time probe Ux p case_1p0 1.0 10 P1 1.234 0.021 case_1p0 1.0 20 P1 1.240 0.019 case_2p0 2.0 10 P1 2.451 0.033 ...为什么用长表而不是宽表长表是pandas/seaborn/ggplot生态的通用语言——一次聚合任意切面出图按参数着色、按监测点分面、按时间变化。如果你一开始就写成每个算例一列的宽表后面每加一个维度都要改数据结构。聚合的三条纪律每一行都要带case标识与参数值——没有来源标识的数据等于没有数据本篇最重要的纪律。时间对齐要显式处理不同算例的时间步可能不同自适应步长聚合时要决定按时间插值对齐还是按步号对齐并把这个决定写在脚本里。缺测要显式标记算例失败或未写出某时刻时用NaN而非静默跳过——静默跳过会让平均值骗人。1.4 pyvista 可视化把场数据变成三维图pyvista是 PyPI 上的三维可视化库基于 VTK 生态适合对场数据做切片、等值面、流线、矢量图。把 OpenFOAM 数据接进pyvista有两条常见路径路径做法优点注意经foamToVTK转换先用 OpenFOAM 工具把场导出为 VTK再用pyvista读稳妥、通用、不依赖 Python 侧解析器会产生额外文件选项以官方文档为准经 Python 读取器用foamlib等直接读场官方 README 展示了加载百万单元的 volVectorField这类能力免中间文件、流程短依赖该库对网格/场的支持范围以官方文档为准推荐姿势先用foamToVTKpyvista路径稳妥、易排错确认流程跑通后再评估是否用 Python 直读以减少中间文件。可视化纪律图必须自解释坐标轴单位、色标范围、物理量名称都要有单位来自算例不能靠记忆色标范围在所有对比图之间保持一致否则看起来更红可能只是色标不同并行/重构后的数据要确认是重构后的第 14 篇——直接对processor*里的场做可视化会得到残缺几何。1.5 出图出表的三条可信任纪律可追溯每个数据点能回溯到哪个算例、哪个文件、哪一行。做法长表必须带case甚至file列。可复现出图脚本与数据在同一份代码库里图表由脚本生成而非手工调整参数如色标范围写进脚本。可审计保留由函数对象产出的原始文件不要只留聚合后的表——原始文件是事实聚合表是观点。二、完整代码与逐行剖析代码 2-1聚合 postProcessing 数据为长表pandas# -*- coding: utf-8 -*- aggregate_postprocessing.py —— 把 postProcessing/ 下的时间序列聚合成对比长表 用法python aggregate_postprocessing.py 算例根目录 [输出csv] 说明先探测目录结构再按文件模式解析不硬编码某个版本的具体路径。 importsysimportrefrompathlibimportPathimportpandasaspddeffind_data_files(case_dir:Path):在 postProcessing/ 下探测数据文件不同函数对象/版本的布局不同故做探测。ppcase_dir/postProcessingifnotpp.is_dir():return[]# 常见为 CSV或数值型数据文件这里按扩展名与文件名模式宽松匹配pats[*.csv,*.dat,*.txt]files[]forpatinpats:files.extend(pp.rglob(pat))returnsorted(files)defparse_time_series(path:Path,case_tag:str)-pd.DataFrame: 解析一个时间序列文件。 注意OpenFOAM 函数对象输出的文件常以 # 开头写表头/元信息 pandas 需要跳过这些注释行comment#。 # 优先按“带表头的表格”解析try:dfpd.read_csv(path,comment#,sepr\s|,,enginepython)exceptException:returnpd.DataFrame()dfdf.loc[:,[cforcindf.columnsifstr(c).strip()!]]ifdf.empty:returndf# 打上来源标识算例、函数对象名父目录、文件路径可追溯纪律df[case]case_tag df[fo_obj]path.parent.parent.name df[source_file]str(path)returndfdefmain()-None:rootPath(sys.argv[1]iflen(sys.argv)1else.)outPath(sys.argv[2]iflen(sys.argv)2elseaggregate.csv)frames[]# 假定 root 下是多个算例目录第 16 篇批量运行产生forcase_dirinsorted(pforpinroot.iterdir()ifp.is_dir()):forfinfind_data_files(case_dir):dfparse_time_series(f,case_tagcase_dir.name)ifnotdf.empty:frames.append(df)ifnotframes:raiseSystemExit(未找到可解析的后处理数据请确认函数对象已产出 postProcessing/ 内容。)all_dfpd.concat(frames,ignore_indexTrue)all_df.to_csv(out,indexFalse)print(f[OK] 已聚合{len(all_df)}行 -{out})print(all_df.head(8).to_string(indexFalse))if__name____main__:main()逐行剖析find_data_files用探测rglob 多模式而非硬编码路径因为函数对象产物的目录层级与文件名随类型/版本不同§一.2。这是不臆造路径的工程做法。pd.read_csv(path, comment#, ...)OpenFOAM 输出文件常以#开头写表头/元信息comment#是 pandas 处理这类文件的标准姿势。sepr\s|,兼容空白或逗号分隔。打上case/fo_obj/source_file三列这是可追溯纪律的直接实现——每行数据都能回溯到具体文件§一.5。解析失败返回空DataFrame而非抛异常批量聚合时个别文件解析失败不应中断全局但也不静默——最终会体现在行数上可加断言。to_csv(out, indexFalse)聚合结果落盘让图可被重生成可复现纪律。代码 2-2多算例对比出图matplotlib风格统一# -*- coding: utf-8 -*- plot_compare.py —— 从聚合长表生成多算例对比图风格统一、自解释 用法python plot_compare.py aggregate.csv importsysimportpandasaspdimportmatplotlib matplotlib.use(Agg)# 无界面环境服务器/CI必须用非交互后端importmatplotlib.pyplotaspltdefmain()-None:csvsys.argv[1]iflen(sys.argv)1elseaggregate.csvdfpd.read_csv(csv)# 1) 基本列检查自解释与可追溯为前提need{case,Time}ifnotneed.issubset(df.columns):raiseSystemExit(f缺少必要列{need}实际列{list(df.columns)})# 2) 选一个数值列作为纵轴按你的函数对象输出列名替换ycols[cforcindf.columnsifcnotin{case,Time,fo_obj,source_file}andpd.api.types.is_numeric_dtype(df[c])]ifnotycols:raiseSystemExit(未找到数值列作为纵轴)ycolycols[0]print(f[INFO] 纵轴列{ycol}如需其它列请按实际列名修改脚本)# 3) 统一风格可复现纪律样式写在脚本里不手工调图plt.rcParams.update({figure.dpi:150,font.size:10,axes.grid:True,grid.alpha:0.3,})fig,axplt.subplots(figsize(7.0,4.2))# 按算例分组画线色标范围在所有对比图之间保持一致本节纪律ymin,ymaxdf[ycol].min(),df[ycol].max()fortag,gindf.groupby(case):gg.sort_values(Time)ax.plot(g[Time],g[ycol],markero,ms3,lw1.2,labeltag)ax.set_xlabel(Time)ax.set_ylabel(ycol)ax.set_ylim(ymin-0.05*abs(ymax-yminor1),ymax0.05*abs(ymax-yminor1))ax.set_title(f{ycol}vs Time多算例对比)ax.legend(titlecase,fontsize8)outcompare.pngfig.tight_layout()fig.savefig(out)print(f[OK] 已生成{out})if__name____main__:main()逐行剖析matplotlib.use(Agg)无界面后端。服务器/CI 上没有显示设备用默认后端会崩——这是自动化出图的必备设置第 20 篇的 CI 会用到。列检查与数值列自动识别让脚本对函数对象列名变化更有鲁棒性同时用[INFO]打印实际使用的列避免悄悄画错列。plt.rcParams统一风格可复现纪律——图表样式是代码的一部分不是事后手调。ax.set_ylim(...)显式设定纵轴范围对应色标/坐标范围在多图间保持一致的纪律这里是线图的等价做法。legend(titlecase) 轴标签带含义图必须自解释§一.4。输出compare.png图由脚本生成任何人重跑都能得到同样的图。代码 2-3场数据可视化foamToVTK pyvista 路径# -*- coding: utf-8 -*- field_slice.py —— 场数据的切片可视化稳妥路径foamToVTK 转换后用 pyvista 读取 用法python field_slice.py 算例目录 时间目录名 前提本机已 source OpenFOAM 环境能执行 foamToVTK并已 pip install pyvista importsubprocessimportsysfrompathlibimportPathdefmain()-None:casePath(sys.argv[1]iflen(sys.argv)1else.)tdirsys.argv[2]iflen(sys.argv)2elselatest# 例如 500# 1) 用 OpenFOAM 工具导出 VTK稳妥、通用具体选项以官方文档为准# 说明-latestTime 取最后时刻具体选项名请用 foamToVTK -help 核对本机版本。cmd[foamToVTK,-case,str(case)]iftdirlatest:cmd.append(-latestTime)else:cmd[-time,tdir]subprocess.run(cmd,checkTrue,cwdcase)print([OK] foamToVTK 完成)# 2) 用 pyvista 读取并做切片importpyvistaaspv vtkssorted((case/VTK).rglob(*.vtk))ifnotvtks:raiseSystemExit(未找到 VTK 文件请确认 foamToVTK 的输出目录以本机输出为准)meshpv.read(vtks[-1])# 取一个代表文件实际项目中应显式选择时刻/区域print(f[INFO] 读入{vtks[-1].name}单元数{mesh.n_cells})# 3) 切片用一个法向平面切过域便于看内部场# 说明origin/normal 需按你的几何设定示例值slmesh.slice(normalz,origin(0.0,0.0,0.005))print(f[INFO] 切片后单元数{sl.n_cells})# 4) 出图无界面后端pv.OFF_SCREENTrueplpv.Plotter(off_screenTrue)pl.add_mesh(sl,show_edgesFalse,cmapviridis)pl.add_axes()pl.camera_positionxy# 俯视 xy 平面按需调整outslice.pngpl.screenshot(out)print(f[OK] 已生成{out}请核对色标与物理量含义后再用于报告)if__name____main__:main()逐行剖析先用foamToVTK导出“稳妥路径”——它不依赖 Python 侧对 OpenFOAM 网格格式的支持范围任何版本都能用。并且我明确提示具体选项名用foamToVTK -help核对本机版本不臆造选项。pv.read(vtks[-1])取代表文件 打印单元数让人一眼看到读进来的到底是什么规模的数据避免读错文件还不知道。mesh.slice(normalz, origin...)切片是看内部场最常用的操作origin我为示例值并标注需按几何设定。pv.OFF_SCREEN Truepl.screenshot(...)无界面出图同Agg后端的意义这是自动化流水线的必要条件。出图后提示请核对色标与物理量含义后再用于报告把可视化纪律写进脚本输出——色标不带单位、含义不明的图是最常见的看起来专业但其实误导的产物。三、常见报错与排查报错 3-1FileNotFoundError/ 聚合结果为空——postProcessing/不存在或为空。现象聚合脚本报未找到可解析的后处理数据。根因算例没配置任何函数对象第 12 篇或函数对象的writeControl/writeInterval导致未写出或算例未运行到写出时刻。解法在controlDict的functions中配置至少一个输出型函数对象如probes运行后确认postProcessing/下有内容也可用postProcess工具对已有时间目录离线补算官方文档说明其无附加选项时执行controlDict中所有时间目录的函数对象。报错 3-2pandas解析出来的列名/数值全是错位的。现象列名是# Time之类或数值列串位。根因OpenFOAM 输出文件的注释行与分隔符与预期不符。解法用comment#跳过注释显式指定分隔符空白或逗号先head看几行原始文件再写解析规则不要凭猜。这是先看数据再写解析的纪律。报错 3-3pyvista读到的几何残缺/只有部分区域。现象可视化结果只有一小块。根因读的是并行processor*目录下的场而不是重构后的串行场第 14 篇。解法先reconstructPar或redistributePar -reconstruct再可视化确认读的路径是算例主目录的时间目录。这是并行算例后处理最常见的坑。报错 3-4多算例对比图看起来差异巨大但数值其实相近。现象图给人和数据不符的印象。根因各图/各曲线的纵轴范围或色标范围不一致。解法在一张图内统一轴范围跨图对比时把色标/轴范围写死在脚本里代码 2-2 的set_ylim在报告里标注轴范围。报错 3-5服务器上出图报 “no display name and no $DISPLAY environment variable”。现象无界面环境出图失败。根因matplotlib/pyvista 尝试打开交互窗口。解法matplotlib 用matplotlib.use(Agg)pyvista 用off_screenTrue并设置pv.OFF_SCREEN True代码 2-2、2-3 已示范。四、动手练习练习 1产出结构化数据在一个官方算例的controlDict中配置probes第 12 篇运行后确认postProcessing/下有文件。判定能指出该文件的具体路径与列含义用head查看原始文件能说明为什么这比解析日志可靠。练习 2长表聚合用代码 2-1 聚合两个以上算例的数据。判定输出 CSV 行数 各算例数据行数之和每行都带case与source_file可追溯能用pandas按case分组统计行数并与原始文件对数。练习 3对比出图用代码 2-2 生成对比图。判定图为多算例曲线、有图例与轴标签、纵轴范围统一图由脚本一键重生成可复现。练习 4场可视化用代码 2-3 对一个已重构的算例做切片出图。判定能看到完整几何不是残缺子域图有色标与坐标轴能说明为什么不能直接可视化processor*下的场。练习 5思考题无标准答案设计一份10 组参数 × 3 个监测点 × 2 个物理量的结果报告方案。验证要点(a) 长表的列设计是否包含case/参数值/监测点/物理量/时间(b) 是否明确时间对齐策略按时间插值还是按步号© 是否处理缺测用 NaN 而非静默跳过(d) 是否保留原始文件以支持审计(e) 出图是否统一轴/色标范围。五、小结与下一篇预告本篇打通了数据到结论的最后一公里三种数据形态postProcessing/时间序列、场数据、日志各有其位——能用函数对象产出结构化数据就不要解析日志postProcess支持离线补算忘记配置也能事后补救pandas长表是聚合的唯一可扩展姿势每行必须带来源标识pyvistafoamToVTK是稳妥的三维可视化路径且必须可视化重构后的数据三条纪律贯穿始终——可追溯、可复现、可审计。第 18 篇《网格与算例流水线自动化》将回到计算之前的环节参数化生成blockMeshDict、snappyHexMesh的分阶段配置、checkMesh质量门禁与失败判定、以及把几何 → 网格 → 质量报告串成全自动流水线——那是把整个前处理过程从手工点选变成一键复现的关键一步。本篇认知问题回显FAQQ1OpenFOAM 的后处理数据有哪几种形态各适合什么工具A三种。第一是时间序列即函数对象输出位于 postProcessing/函数对象名/时间/ 下的文件常为 CSV结构化且按时间排列最适合 pandas 聚合、趋势对比与出图。第二是场数据位于各时间目录下如 500/U、500/p为 OpenFOAM 自有文本或二进制场格式最适合三维可视化与空间分析。第三是日志文本log.求解器非结构化只适合快速排查与残差轨迹不适合作为长期数据源。关键判断是能用函数对象产出结构化数据就不要解析日志因为日志格式随版本与求解器演进基于日志的解析脚本最脆弱。Q2postProcessing/ 下的文件是什么结构为什么比解析日志可靠A函数对象产物一般落在 postProcessing/ 下按函数对象名、时间、文件名分层组织但其目录层级与文件命名随函数对象类型与版本不同因此解析脚本应先探测目录结构再按模式匹配文件不要硬编码路径。它比日志可靠的原因是结构固定、由官方函数对象写入、含明确的时间列与物理量列而日志是非结构化文本且随版本与求解器变化。此外官方 postProcess 工具在没有附加选项时会执行 controlDict 中为所有时间目录列出的函数对象意味着可对已有时间目录离线补算算例跑完后发现忘了加监测也不必重跑。Q3怎么用 pandas 把多算例、多时刻的数据聚合成对比表A用长表long format列为 case、参数值、time、监测点、物理量等每个数据点一行。做法是用 pandas 读各算例 postProcessing/ 下的文件用 comment‘#’ 跳过 OpenFOAM 的注释行显式指定分隔符再打上来源标识列case、函数对象名、源文件路径后 concat。之所以用长表而非宽表是因为长表是 pandas 与绘图生态的通用语言一次聚合即可任意切面出图而宽表每加一个维度都要改数据结构。三条纪律每行必须带来源标识、时间对齐策略要显式决定按时间插值或按步号、缺测用 NaN 显式标记而非静默跳过。Q4怎么把场数据读进 pyvista 做三维可视化A两条路径。稳妥路径是先用 OpenFOAM 的 foamToVTK 把场导出为 VTK再用 pyvista 读取不依赖 Python 侧对网格格式的支持范围易排错具体命令行选项以官方文档与本机 help 为准。短流程路径是用 foamlib 等库直接读场官方 README 展示了加载百万单元 volVectorField 这类能力但依赖库对网格与场的支持范围。推荐先跑通 foamToVTK 加 pyvista再评估是否减少中间文件。关键纪律是可视化必须使用重构后的数据直接对 processor* 目录下的场出图会得到残缺几何因此应先 reconstructPar 或 redistributePar -reconstruct。Q5出图出表有哪些可信任纪律A三条。可追溯每个数据点能回溯到哪个算例、哪个文件、哪一行因此长表必须带 case 与来源文件列。可复现出图脚本与数据同在一份代码库图表由脚本生成而非手工调整坐标轴范围与色标范围等参数写进脚本无界面环境需使用 matplotlib 的 Agg 后端或 pyvista 的 off_screen 模式。可审计保留函数对象产出的原始文件不要只留聚合后的表因为原始文件是事实、聚合表是观点。此外图必须自解释轴标签、单位、色标含义齐全多图对比时色标与轴范围应保持一致否则会出现看起来差异巨大而数值其实相近的误导。
返回列表