How do I calculate the numerical difference between two fields stored in two different VTK files with the same structure?

Suppose I have two VTK files, both in structured grid format. The structured grids are the same (they have the same list of points, in the same order), and there is a field, call it "Phi", in each VTK file. I want to create a third VTK file, again with the same structured grid, and plot a field that is the difference between Phi in the first VTK file and Phi in the second VTK file.

I know how to do this manually; I can parse the raw text in the two VTK files, copy the data into arrays, subtract one array from the other, and then dump the data in the right format into a new file. Is there a better way to calculate this difference and export it to VTK? A solution in Python, or in visualization software like VisIt or Paraview would be preferable to using a compiled language like C++.

The purpose of calculating this difference is to compare different numerical methods for calculating the solution of a PDE; since I'm using the same software to generate the solutions, I can guarantee that all data but the field Phi will be the same in each file I generate.

• I posted this question because it took me around a day and a half to figure out the answer; had I not found it yesterday, I would've asked this question on here anyway. I am interested in seeing if there are other quick ways of accomplishing the same task. Nov 7 '13 at 19:01
• When you say "parse the raw text" you mean literally go into the file or use a python parser ?
Sep 2 '14 at 13:09
• At the time, I meant write a Python parser by hand. Sep 2 '14 at 14:21

The easiest way I could find to subtract two fields from different VTK files with the same structured grid is to use a programmable filter in Paraview, which lets you manipulate data using Python scripts.

In the programmable filter dialogue box, you can subtract the two arrays and write to output with the code:

   phi_0 = inputs[0].CellData['Phi']
phi_1 = inputs[1].CellData['Phi']
output.CellData.append(phi_1 - phi_0, 'difference')


In this instance, the field Phi happens to be cell data. If your field is point data, replace CellData everywhere in the script with PointData. See http://public.kitware.com/pipermail/paraview/2010-April/016667.html for further details.

• It is never too much to remind that to have two inputs (inputs[0] and inputs[1]), one must highlight both data sets before selecting the Programmable Filter (this is mentioned one the referred link). Jan 18 '16 at 14:29

In ParaView there is the Append Attributes filter which can be used for this. It requires that the same number of points are in the data set for appending point data properly and the same number of cells are in the data set for appending cell data properly. It will have problems though with arrays of the same name (i.e. Phi in your example). You can copy that array easily with the Calculator filter though before using the Append Attributes filter. Then you can use another Calculator filter to do the subtraction. This is probably less efficient than using ParaView's Python programmable filter though. Besides that, you could use the vtkpython executable to do it manually since you'd have direct access to both grids and their attributes.

I don't have a particularly good approach, but I'd copy the 'phi' field from one VTK file into the other and name it 'phiprime' or something. In both Paraview and Visit, you have the option of defining new fields by a formula that uses the values of other fields. You can then define a field "error" as "error=phi-phiprime" in the field editor, and plot this field "error" as either a surface, a contour plot, or whatever you're interested in.

The step of copying the block of data from one file to the other is clearly awkward, but it's the best I can come up with.

I realized this is a bit older, but I though you might be interested in the VisIt solution:

You can do this in VisIt with something called a Connectively Based Cross Mesh Field Expression. That is a mouthful, but is basically machinery to maps fields between databases (in your case VTK files).

The "Connectivity based" (conn_cmfe) is used when the topology is the same between the files -- like in your case.

There is also a "Position based" (pos_cmfe) that samples between meshes with different topologies.

For your case, open the first file & use the Expressions window to define an expression (MyPhi_Diff) :

Phi - conn_cmfe(<file2.vtk:Phi>, mesh)


Then you can plot "MyPhi_Diff" with a Pseudocolor plot.

There is also wizard that you can use to help define the expression