前面的博客中已经提到，Python 有 Matplotlib
这种强大的包可以「在线」地生成(并保存)漂亮的图形，但对于计算流体力学的后处理来说，很多时候需要更多「事后」的后处理。在线的即时输出，虽然可以快速预览结果，但对于想输出什么结果必须在写程序阶段就全部想明白，而不是先计算完，按一定间隔输出整个流场的信息，然后后期想要什么数据慢慢提取分析就行。
解决这个问题的方法也很简单，直接照抄开源 LBM 软件 Palabos 的方案，用 vtk
格式保存流场信息，然后用 Paraview 进行后处理。
这个「十分简单」的技术路线结果却花了我很多天的时间，才基本搞清楚其中的奥秘，从多个 tricky
part 中脱身。技术上来说，Python 输出 vtk 有成熟的软件包(我采用的是pyevtk
这个库[1])，实现方法也比较直观。例如官方提供的一个示例是这样的(我稍作了修改)：
from evtk.vtk import VtkFile,
VtkRectilinearGrid
import numpy as np
nx, ny, nz = 6, 6, 2
lx, ly, lz = 1.0, 1.0, 1.0
dx, dy, dz = lx/nx, ly/ny, lz/nz
npoints = (nx + 1) * (ny + 1) * (nz +
1)
x = np.arange(0, lx + 0.1*dx, dx,
dtype='float64')
y = np.arange(0, ly + 0.1*dy, dy,
dtype='float64')
z = np.arange(0, lz + 0.1*dz, dz,
dtype='float64')
start, end = (0,0,0), (nx, ny, nz)
w = VtkFile("./evtk_test",
VtkRectilinearGrid)
w.openGrid(start = start, end =
end)
w.openPiece( start = start, end =
end)
# Point data
temp = np.random.rand(npoints)
vx = vy = vz = np.zeros([nx + 1, ny + 1, nz +
1], dtype="float64", order = 'F')
w.openData("Point", scalars = "Temperature",
vectors = "Velocity")
w.addData("Temperature", temp)
w.addData("Velocity", (vx,vy,vz))
w.closeData("Point")
# Cell data
# 因为 LBM 计算中结果都在 point 上，这里省略。
# Coordinates of cell vertices
w.openElement("Coordinates")
w.addData("x_coordinates", x);
w.addData("y_coordinates", y);
w.addData("z_coordinates", z);
w.closeElement("Coordinates");
w.closePiece()
w.closeGrid()
w.appendData(data = temp)
w.appendData(data = (vx,vy,vz))
w.appendData(x).appendData(y).appendData(z)
w.save()
这里真正关于 vtk 输出的逻辑非常简单，基本上就跟 Python
写入普通文件的流程一样，唯一不同的是这里似乎有点冗余，先要 addData 一下，然后又要 appendData。
但是，在自己写程序的时候，却碰到了一些问题。
首先的问题是，这个示例中展示的是3维的问题，而我初期 LBM 的计算都是2维的，怎么办？
很容易想到的办法是，把2维的数据映射到 z=0 平面上，然后再复制几份到 z=1,2,3…nz
平面上，这样输出一个沿z轴完全相同的伪3维数据。办法虽然土了点，但简单有效！当然，后来发现更简单且不这么丑陋的办法是，只输出到 z=0
的平面上即可，即完全按3维的情况输出，只是z方向只有一组数据而已。这样导入 Paraview 时，可以直接就被识别成2维结果。
Vtk 文件输出了，进入 Paraview
也能正常打开，显示温场，正常，显示流场，全乱了！怎么回事？
经过至少3天断断续续的排查，总算弄明白了其中的蹊跷，全在一个非常小的地方：附加的数据不能用切片。我的速度数据是存储在一个 u =
np.array([nx+1, ny+1, 1, 3]) 的数组中，于是在 appendData
中就把三个速度分量写成了切片的形式：
ux, uy, uz = u_3d[:,:,:,0], u_3d[:,:,:,1], u_3d[:,:,:,2]
结果这样输出的结果就怎么都不对。最后发现，只要不再用切片输出数据，而是复制一下上面的切片(np.copy
方法)，结果就正常了。
实在是莫名其妙的问题，更莫名其妙的解决方案。
倒过头来看，这个问题肯定是由于存储数据的顺序决定的(后面是一段不太靠谱的分析)。u、t
等多维数组，虽然本身是多维的，但在存储的时候，肯定还是展开成一维的形式存储的。而多维数据的存储，有C和Fortran两种不同的顺序(注意
pyevtk 给出的示例中用的是 order='F')。Pyevtk
库可能内置了什么判断，可以知道输入的数组是什么类型的排序方式，但如果输入的是切片，就扰乱了这个判断，从而导致写入的顺序错误。
PS. 补充一下，之前用的复制多个z数据，强行输出伪3维结果的方法中，如果只输出两层结果，在
Paraview 中，可以正常显示标量云图和矢量图，但无法显示流线图，必须复制多份(我自己的测试是4份就行了)。
[1] Python 似乎还有专门的库用于处理 vtk 格式，如 tvtk 等