Fluidity

Merge lp:~fluidity-core/fluidity/speedup-vtudiff into lp:fluidity

  1. speedup-vtudiff
  2. Merge into dev-trunk
Proposed by Stephan Kramer
Status: Merged
Merged at revision: 4320
Proposed branch: lp:~fluidity-core/fluidity/speedup-vtudiff
Merge into: lp:fluidity
Diff against target: 386 lines (+101/-212)
2 files modified
manual/visualisation_and_diagnostics.tex (+1/-13)
python/vtktools.py (+100/-199)
To merge this branch: bzr merge lp:~fluidity-core/fluidity/speedup-vtudiff
Reviewer Review Type Date Requested Status
Jon Hill Approve
Review via email: mp+201957@code.launchpad.net

Description of the change

This speeds up the vtudiff script run on vtus with multiple fields, by only setting up the VTK probe once instead of once per field. It does not change anything in the way fields are probed.

There is also a little bit of cleaning up:
* removing the option to call vtktools as a main: this really doesn't make sense and should be in a separate script
* removing the "field manipulation" methods. These were basically vectorial operations on fields but implemented in a rather inefficient way. I see no reason to not simply use numpy vector operations for that, so instead of:

vtu.SubFieldFromField("fieldName", array, "newFieldName")

just do:

field = vtu.GetField("fieldName")
vtu.AddField("newFieldName", field-array)

To post a comment you must log in.
lp:~fluidity-core/fluidity/speedup-vtudiff updated
4308. By Stephan Kramer

For vtudiff: Also compute diff of cell-based fields (only works if both meshes are the same)

Revision history for this message
Jon Hill (jon-hill) :
review: Approve
Revision history for this message
Stephan Kramer (s-kramer) wrote :

Thanks Jon, will merge once I have a green buildbot.

lp:~fluidity-core/fluidity/speedup-vtudiff updated
4309. By Stephan Kramer

Maintain the same sign convention in vtudiff as before:

With "vtudiff A.vtu B.vtu C.vtu" we get fieldC = fieldA-fieldB
Also don't try to subtract vtkGhostLevels field.

4310. By Stephan Kramer

Merge in trunk.

Preview Diff

[H/L] Next/Prev Comment, [J/K] Next/Prev File, [N/P] Next/Prev Hunk
1=== modified file 'manual/visualisation_and_diagnostics.tex'
2--- manual/visualisation_and_diagnostics.tex 2013-04-12 10:31:41 +0000
3+++ manual/visualisation_and_diagnostics.tex 2014-02-19 10:25:38 +0000
4@@ -350,8 +350,6 @@
5 method & arguments & use \\ \hline
6 \lstinline[language=Python]+AddField+ & \lstinline[language=Python]+name, array+ & Adds the values in \lstinline[language=Python]+array+ (the entries of which may have an arbitrary number of components) as a field called \lstinline[language=Python]+name+ \\ \hline
7 %
8-\lstinline[language=Python]+AddFieldtoField+ & \lstinline[language=Python]+fieldname,+ \lstinline[language=Python]+array,+ \lstinline[language=Python]+newFieldName=None+ & Adds the values in \lstinline[language=Python]+array+ to the field \lstinline[language=Python]+fieldname+. If \lstinline[language=Python]+newFieldName+ is specified then a new field with that name is created with the new values, otherwise the original field is replaced. \\ \hline
9-%
10 \lstinline[language=Python]+AddScalarField+ & \lstinline[language=Python]+name, array+ & Adds a scalar field called \lstinline[language=Python]+name+ using the values in \lstinline[language=Python]+array+ \\ \hline
11 %
12 \lstinline[language=Python]+AddVectorField+ & \lstinline[language=Python]+name, array+ & Adds a vector field called \lstinline[language=Python]+name+ using the values in \lstinline[language=Python]+array+ \\ \hline
13@@ -365,10 +363,6 @@
14 \lstinline[language=Python]+CellDataToPointData+ & & transforms all cell--wise fields to point--wise fields. All existing fields will remain, \\ \hline
15 \lstinline[language=Python]+Crop+ & \lstinline[language=Python]+min_x, max_x,+ \lstinline[language=Python]+min_y, max_y+, \lstinline[language=Python]+min_z, max_z+ & Crops the edges defined by the bounding box given by the arguments \\ \hline
16 %
17-\lstinline[language=Python]+CrossFieldWithField+ & \lstinline[language=Python]+fieldName,+ \lstinline[language=Python]+array,+ \lstinline[language=Python]+newFieldName=None,+ \lstinline[language=Python]+postMultiply=True+ & Calculates the cross product, ${\bf a} \times {\bf b}$, where ${\bf a}$ is the field \lstinline[language=Python]+fieldName+ and ${\bf b}$ is \lstinline[language=Python]+array+. If \lstinline[language=Python]+postMultiply$\neq$True+ then ${\bf b} \times {\bf a}$ will be calculated. If \lstinline[language=Python]+newFieldName+ is specified then a new field with that name is created that takes the values of the cross product, otherwise the original field is replaced. \\ \hline
18-%
19-\lstinline[language=Python]+DotFieldWithField+ & \lstinline[language=Python]+fieldName,+ \lstinline[language=Python]+array,+ \lstinline[language=Python]+newFieldName+ & Calculates the dot product of the field called \lstinline[language=Python]+fieldName+ and \lstinline[language=Python]+array+. If \lstinline[language=Python]+newFieldName+ is specified then a new field with that name is created that takes the values of the dot product, otherwise the original field is replaced. \\ \hline
20-%
21 \lstinline[language=Python]+GetCellPoints+ & \lstinline[language=Python]+id+ & Returns an array with the node numbers of the cell (mesh element) number \lstinline[language=Python]+id+. \\ \hline
22 %
23 \lstinline[language=Python]+GetCellVolume+ & \lstinline[language=Python]+id+ & Returns the volume of the cell (mesh element) with number \lstinline[language=Python]+id+ \\ \hline
24@@ -404,18 +398,12 @@
25 %
26 \lstinline[language=Python]+IntegrateField+ & \lstinline[language=Python]+field+ & Returns the integral of the field called \lstinline[language=Python]+field+ assuming a linear representation on a tetrahedral mesh. \\ \hline
27 %
28-\lstinline[language=Python]+ManipulateField+ & \lstinline[language=Python]+fieldName,+ \lstinline[language=Python]+manipFunc,+ \lstinline[language=Python]+newFieldName=None+ & Generic field manipulation method. Applies the supplied manipulation function, \lstinline[language=Python]+manipFunc+, to the field called \lstinline[language=Python]+fieldName+. \lstinline[language=Python]+manipFunc+ must have form \lstinline[language=Python]+def manipFunc(field, index):+ $\ldots$ \lstinline[language=Python]+return fieldValAtIndex+. If \lstinline[language=Python]+newFieldName+ is specified then a new field with that name is created that takes the calculated values, otherwise the original field is replaced. \\ \hline
29-%
30-\lstinline[language=Python]+MatMulFieldWithField+ & \lstinline[language=Python]+fieldName,+ \lstinline[language=Python]+array,+ \lstinline[language=Python]+newFieldName=None,+ \lstinline[language=Python]+postMultiply=True+. & Multiplies two matrices $\bar{\bar{A}}\bar{\bar{B}}$ where $\bar{\bar{A}}$ is the field \lstinline[language=Python]+fieldName+, and $\bar{\bar{B}}$ is \lstinline[language=Python]+array+. If \lstinline[language=Python]+postMultiply+$\neq$\lstinline[language=Python]+True+ then $\bar{\bar{A}}\bar{\bar{B}}$ will be calculated. If \lstinline[language=Python]+newFieldName+ is specified then a new field with that name is created that takes the values of the product, otherwise the original field is replaced. \\ \hline
31-%
32- \lstinline[language=Python]+ProbeData+ & \lstinline[language=Python]+coordinates,+ \lstinline[language=Python]+name+ & Returns an array of values of the field called \lstinline[language=Python]+name+ at the positions given in \lstinline[language=Python]+coordinates+. The values are calculated by interpolation of the field to the positions given. \lstinline[language=Python]+coordinates+ can be created using \lstinline[language=Python]+vtktools.arr()+ e.g. \lstinline[language=Python]+coordinates=vtktools.arr([[1,1,1],[1,1,2]])+. \\ \hline
33+\lstinline[language=Python]+ProbeData+ & \lstinline[language=Python]+coordinates,+ \lstinline[language=Python]+name+ & Returns an array of values of the field called \lstinline[language=Python]+name+ at the positions given in \lstinline[language=Python]+coordinates+. The values are calculated by interpolation of the field to the positions given. \lstinline[language=Python]+coordinates+ can be created using \lstinline[language=Python]+vtktools.arr()+ e.g. \lstinline[language=Python]+coordinates=vtktools.arr([[1,1,1],[1,1,2]])+. \\ \hline
34 %
35 \lstinline[language=Python]+RemoveField+ & \lstinline[language=Python]+name+ & Removes the field called \lstinline[language=Python]+name+. \\ \hline
36 %
37 \lstinline[language=Python]+StructuredPointProbe+ & \lstinline[language=Python]+nx,+ \lstinline[language=Python]+ny,+ \lstinline[language=Python]+nz,+ \lstinline[language=Python]+bounding_box=None+ & Returns a vtk structured object. \lstinline[language=Python]+nx, ny, nz+ are the number of points in the $x,\, y, \, z$ directions respectively. If \lstinline[language=Python]+bounding_box+ is not specified the bounding box of the domain is calculated automatically. If specified \lstinline[language=Python]+bounding_box = [xmin, xmax, ymin, ymax, zmin, zmax]+. \\ \hline
38 %
39-\lstinline[language=Python]+SubFieldFromField+ & \lstinline[language=Python]+fieldName,+ \lstinline[language=Python]+array,+ \lstinline[language=Python]+newFieldName=None+ & Subtracts \lstinline[language=Python]+array+ from the field called \lstinline[language=Python]+fieldName+. If \lstinline[language=Python]+newFieldName+ is specified then a new field with that name is created that takes the values of the product, otherwise the original field is replaced. \\ \hline
40-%
41 \lstinline[language=Python]+Write+ & \lstinline[language=Python]+filename = []+ & Writes the data to a vtu file. If \lstinline[language=Python]+filename+ is not specified the name of the file originally read in will be used and therefore the input file will be overwritten. \\ \hline
42 %
43 %
44
45=== modified file 'python/vtktools.py'
46--- python/vtktools.py 2012-03-19 21:00:50 +0000
47+++ python/vtktools.py 2014-02-19 10:25:38 +0000
48@@ -262,64 +262,8 @@
49
50 def ProbeData(self, coordinates, name):
51 """Interpolate field values at these coordinates."""
52-
53- # Initialise locator
54- locator = vtk.vtkPointLocator()
55- locator.SetDataSet(self.ugrid)
56- locator.SetTolerance(10.0)
57- locator.Update()
58-
59- # Initialise probe
60- points = vtk.vtkPoints()
61- points.SetDataTypeToDouble()
62- ilen, jlen = coordinates.shape
63- for i in range(ilen):
64- points.InsertNextPoint(coordinates[i][0], coordinates[i][1], coordinates[i][2])
65- polydata = vtk.vtkPolyData()
66- polydata.SetPoints(points)
67- probe = vtk.vtkProbeFilter()
68- probe.SetInput(polydata)
69- probe.SetSource(self.ugrid)
70- probe.Update()
71-
72- # Generate a list invalidNodes, containing a map from invalid nodes in the
73- # result to their closest nodes in the input
74- valid_ids = probe.GetValidPoints()
75- valid_loc = 0
76- invalidNodes = []
77- for i in range(ilen):
78- if valid_ids.GetTuple1(valid_loc) == i:
79- valid_loc += 1
80- else:
81- nearest = locator.FindClosestPoint([coordinates[i][0], coordinates[i][1], coordinates[i][2]])
82- invalidNodes.append((i, nearest))
83-
84- # Get final updated values
85- pointdata=probe.GetOutput().GetPointData()
86- vtkdata=pointdata.GetArray(name)
87- nc=vtkdata.GetNumberOfComponents()
88- nt=vtkdata.GetNumberOfTuples()
89- array = arr([vtkdata.GetValue(i) for i in range(nt * nc)])
90-
91- # Fix the point data at invalid nodes
92- if len(invalidNodes) > 0:
93- try:
94- oldField = self.ugrid.GetPointData().GetArray(name)
95- components = oldField.GetNumberOfComponents()
96- except:
97- try:
98- oldField = self.ugrid.GetCellData().GetArray(name)
99- components = oldField.GetNumberOfComponents()
100- except:
101- raise Exception("ERROR: couldn't find point or cell field data with name "+name+" in file "+self.filename+".")
102- for invalidNode, nearest in invalidNodes:
103- for comp in range(nc):
104- array[invalidNode * nc + comp] = oldField.GetValue(nearest * nc + comp)
105-
106- valShape = self.GetField(name)[0].shape
107- array.shape = tuple([nt] + list(valShape))
108-
109- return array
110+ probe = VTU_Probe(self.ugrid, coordinates)
111+ return probe.GetField(name)
112
113 def RemoveField(self, name):
114 """Removes said field from the unstructured grid."""
115@@ -489,93 +433,6 @@
116 probe.Update ()
117
118 return probe.GetOutput()
119-
120- ### Field manipulation methods ###
121-
122- def ManipulateField(self, fieldName, manipFunc, newFieldName = None):
123- """
124- Generic field manipulation method. Applies the supplied manipulation function
125- manipFunc to the field fieldName. manipFunc must be a function of the form:
126-
127- def manipFunc(field, index):
128- # ...
129- return fieldValAtIndex
130- """
131-
132- field = self.GetField(fieldName)
133- if newFieldName is None or fieldName == newFieldName:
134- self.RemoveField(fieldName)
135- newFieldName = fieldName
136-
137- field = arr([manipFunc(field, i) for i in range(len(field))])
138- self.AddField(newFieldName, field)
139-
140- return
141-
142- def AddFieldToField(self, fieldName, array, newFieldName = None):
143- def ManipFunc(field, index):
144- return field[index] + array[index]
145-
146- self.ManipulateField(fieldName, ManipFunc, newFieldName)
147-
148- return
149-
150- def SubFieldFromField(self, fieldName, array, newFieldName = None):
151- def ManipFunc(field, index):
152- return field[index] - array[index]
153-
154- self.ManipulateField(fieldName, ManipFunc, newFieldName)
155-
156- return
157-
158- def DotFieldWithField(self, fieldName, array, newFieldName = None):
159- """
160- Dot product
161- """
162-
163- def ManipFunc(field, index):
164- sum = 0.0
165- for i, val in enumerate(field[index]):
166- sum += val * array[index][i]
167-
168- return sum
169-
170- self.ManipulateField(fieldName, ManipFunc, newFieldName)
171-
172- return
173-
174- def CrossFieldWithField(self, fieldName, array, newFieldName = None, postMultiply = True):
175- """
176- Cross product
177- """
178-
179- def ManipFunc(field, index):
180- if postMultiply:
181- return numpy.cross(field[index], array[index])
182- else:
183- return numpy.cross(array[index], field[index])
184-
185- self.ManipulateField(fieldName, ManipFunc, newFieldName)
186-
187- return
188-
189- def MatMulFieldWithField(self, fieldName, array, newFieldName = None, postMultiply = True):
190- """
191- Matrix multiplication
192- """
193-
194- def ManipFunc(field, index):
195- if postMultiply:
196- return numpy.matrix(field[i]) * numpy.matix(array[i])
197- else:
198- return numpy.matix(array[i]) * numpy.matrix(field[i])
199-
200- self.ManipulateField(fieldName, ManipFunc, newFieldName)
201-
202- return
203-
204- # Default multiplication is dot product
205- MulFieldByField = DotFieldWithField
206
207 def GetDerivative(self, name):
208 """
209@@ -629,6 +486,73 @@
210 cdtpd.PassCellDataOn()
211 cdtpd.Update()
212 self.ugrid=cdtpd.GetUnstructuredGridOutput()
213+
214+class VTU_Probe(object):
215+ """A class that combines a vtkProbeFilter with a list of invalid points (points that it failed to probe
216+ where we take the value of the nearest point)"""
217+
218+ def __init__(self, ugrid, coordinates):
219+ # Initialise locator
220+ locator = vtk.vtkPointLocator()
221+ locator.SetDataSet(ugrid)
222+ locator.SetTolerance(10.0)
223+ locator.Update()
224+
225+ # Initialise probe
226+ points = vtk.vtkPoints()
227+ points.SetDataTypeToDouble()
228+ ilen, jlen = coordinates.shape
229+ for i in range(ilen):
230+ points.InsertNextPoint(coordinates[i][0], coordinates[i][1], coordinates[i][2])
231+ polydata = vtk.vtkPolyData()
232+ polydata.SetPoints(points)
233+ self.probe = vtk.vtkProbeFilter()
234+ self.probe.SetInput(polydata)
235+ self.probe.SetSource(ugrid)
236+ self.probe.Update()
237+
238+ # Generate a list invalidNodes, containing a map from invalid nodes in the
239+ # result to their closest nodes in the input
240+ valid_ids = self.probe.GetValidPoints()
241+ valid_loc = 0
242+ self.invalidNodes = []
243+ for i in range(ilen):
244+ if valid_ids.GetTuple1(valid_loc) == i:
245+ valid_loc += 1
246+ else:
247+ nearest = locator.FindClosestPoint([coordinates[i][0], coordinates[i][1], coordinates[i][2]])
248+ self.invalidNodes.append((i, nearest))
249+ self.ugrid = ugrid
250+
251+ def GetField(self, name):
252+ # Get final updated values
253+ pointdata = self.probe.GetOutput().GetPointData()
254+ vtkdata=pointdata.GetArray(name)
255+ nc=vtkdata.GetNumberOfComponents()
256+ nt=vtkdata.GetNumberOfTuples()
257+ array = arr([vtkdata.GetValue(i) for i in range(nt * nc)])
258+
259+ # Fix the point data at invalid nodes
260+ if len(self.invalidNodes) > 0:
261+ oldField = self.ugrid.GetPointData().GetArray(name)
262+ if oldField is None:
263+ oldField = self.ugrid.GetCellData().GetArray(name)
264+ if oldField is None:
265+ raise Exception("ERROR: couldn't find point or cell field data with name "+name+".")
266+ components = oldField.GetNumberOfComponents()
267+ for invalidNode, nearest in self.invalidNodes:
268+ for comp in range(nc):
269+ array[invalidNode * nc + comp] = oldField.GetValue(nearest * nc + comp)
270+
271+ # this is a copy and paster from vtu.GetField above:
272+ if nc==9:
273+ return array.reshape(nt,3,3)
274+ elif nc==4:
275+ return array.reshape(nt,2,2)
276+ else:
277+ return array.reshape(nt,nc)
278+
279+ return array
280
281 def VtuMatchLocations(vtu1, vtu2, tolerance = 1.0e-6):
282 """
283@@ -698,6 +622,8 @@
284
285 # If the input vtu point locations match, do not use probe
286 useProbe = not VtuMatchLocations(vtu1, vtu2)
287+ if useProbe:
288+ probe = VTU_Probe(vtu2.ugrid, vtu1.GetLocations())
289
290 # Copy the grid from the first input vtu into the output vtu
291 resultVtu.ugrid.DeepCopy(vtu1.ugrid)
292@@ -707,65 +633,40 @@
293 fieldNames1 = vtu1.GetFieldNames()
294 fieldNames2 = vtu2.GetFieldNames()
295 for fieldName in fieldNames1:
296+ field1 = vtu1.GetField(fieldName)
297 if fieldName in fieldNames2:
298 if useProbe:
299- field2 = vtu2.ProbeData(vtu1.GetLocations(), fieldName)
300+ field2 = probe.GetField(fieldName)
301 else:
302 field2 = vtu2.GetField(fieldName)
303- resultVtu.SubFieldFromField(fieldName, field2)
304+ resultVtu.AddField(fieldName, field1-field2)
305 else:
306 resultVtu.RemoveField(fieldName)
307
308+ # Also look for cell-based fields. This only works if we don't have
309+ # to interpolate (both meshes are the same)
310+ vtkdata=vtu1.ugrid.GetCellData()
311+ fieldNames1 = [vtkdata.GetArrayName(i) for i in range(vtkdata.GetNumberOfArrays())]
312+ vtkdata=vtu2.ugrid.GetCellData()
313+ fieldNames2 = [vtkdata.GetArrayName(i) for i in range(vtkdata.GetNumberOfArrays())]
314+ if useProbe:
315+ # meshes are different - we can't interpolate cell-based fields so let's just remove them from the output
316+ for fieldName in fieldNames1:
317+ if fieldName=='vtkGhostLevels':
318+ # this field should just be passed on unchanged
319+ continue
320+ resultVtu.RemoveField(fieldName)
321+ else:
322+ # meshes are the same - we can simply subtract
323+ for fieldName in fieldNames1:
324+ if fieldName=='vtkGhostLevels':
325+ # this field should just be passed on unchanged
326+ continue
327+ elif fieldName in fieldNames2:
328+ field1 = vtu1.GetField(fieldName)
329+ field2 = vtu2.GetField(fieldName)
330+ resultVtu.AddField(fieldName, field1-field2)
331+ else:
332+ resultVtu.RemoveField(fieldName)
333+
334 return resultVtu
335-
336-def usage():
337- print 'Usage:'
338- print 'COMMAND LINE: vtktools.py [-h] [-p] [-e var1,var2, ...] INPUT_FILENAME'
339- print ''
340- print 'INPUT_FILENAME:'
341- print ' The input file name.'
342- print ''
343- print 'OPTIONS:'
344- print ' -h Prints this usage message.'
345- print ' -p Converts the coordinates from xyz to latitude and longitude.'
346- print ' -e Extracts the data point from the variables provided.'
347-
348-if __name__ == "__main__":
349- import vtktools
350- import math
351- import getopt
352-
353- optlist, args = getopt.getopt(sys.argv[1:], 'hpe:')
354-
355- v = vtktools.vtu(args[0])
356-
357- # Parse arguments
358- LongLat = False
359- for o,a in optlist:
360- if o == '-h':
361- usage()
362- elif o == '-p':
363- LongLat = True
364- elif o == '-e':
365- scalars = a.strip().split(",")
366-
367- # Project domain if necessary
368- if(LongLat):
369- v.ApplyEarthProjection()
370-
371- # Extract variables
372- if(scalars):
373- npoints = v.ugrid.GetNumberOfPoints ()
374- nvar = len(scalars)
375- for i in range (npoints):
376- (x,y,z) = v.ugrid.GetPoint (i)
377- line = "%lf "%x+"%lf "%y
378- for scalar in scalars:
379- line = line+" %lf"%v.ugrid.GetPointData().GetArray(scalar).GetTuple1(i)
380- print line
381-
382-
383-
384-
385-
386-