In this tutorial you will learn how to use turbulucid to analyse and plot data. We will use a dataset from a computational fluid dynamics simulation of a flow over a backward-facing step (BFS) as an example to learn about the different functionalities of the package!
Lets start by doing some imports and making figures a bit larger than the default size.
In [1]:
%matplotlib inline
In [2]:
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from turbulucid import *
matplotlib.rcParams['figure.figsize'] = (12, 6)
The most important class defined by turbulucid is Case. A Case object contains the dataset that is being analysed and several conenience function for data extraction. Opening a dataset therefore amounts to constructing a Case object, which requires a path to the file with the data.
The file can be either VTK polyData in the legacy format (.vtk), VTK polyData in XML format (.vtu) or a multiblock VTK dataset (.vtm) where one block corresponds to the internal data and the other blocks are boundaries.
Let us open the BFS dataset stored as a .vtm.
In [3]:
from os.path import join
import turbulucid
case = Case(join(turbulucid.__path__[0], "datasets", "bfs", "bfs.vtm"))
The dataset is now opened. Let us inspect some of its properties. The fields attribute returns a list of names of the CellData arrays that exist in the dataset.
In [4]:
print(case.fields)
We can access each field by using the Case object as a dictionary, i.e Case['UMean'] will return the values of UMean as a numpy array.
In [5]:
case['UMean']
Out[5]:
The order of the data corresponds to the order of the cells in the dataset. This doesn't necessarily follow any logical pattern.
However, it is possible to get the location of each cell centre via Case.cellCentres.
In [6]:
case.cellCentres
Out[6]:
The order in the cellCentres array is the same as that of the data for each field. This means that it is possible, for intance, to find the value of field at a certian location by first finding the closest cell centre and the grabbing the field's value at the same index.
It is possible to add additional fields to the case. Let us create a new field -- the magnitude of the mean velocity UMean.
In [7]:
case["magUMean"] = np.linalg.norm(case["UMean"], axis=1)
case["magUMean"]
Out[7]:
It is time to make a plot! Let us plot the newly created field. This can be done by using the plot_field function. The first argument to the function is a Case object and the second is the field we want to plot.
In [8]:
plot_field(case, "magUMean")
Out[8]:
Voila! Note that the returned object of plot_field is a PatchCollection. This means that by assigning plot_field to a variable you can further customize the figure. Also, you can pass additional keyword arguments to plot_field that will be used in the PatchCollection constructor. This allows for a lot of customization. Let us, for instance, use a different colormap.
In [9]:
plot_field(case, 'magUMean', cmap='Blues')
Out[9]:
Sometimes it is nice to scale the axes by a certain parameter. For the BFS it is common to scale the axes with the height of the step h. To this end plot_field and all the other plot functions in turbulucid accept the keyword arguments scaleX and scaleY.
In [10]:
h = 0.01
plot_field(case, 'magUMean', scaleX=h, scaleY=h, cmap='Blues')
plt.xlabel(r'$x/h$')
plt.ylabel(r'$y/h$')
Out[10]:
Don't forget that all function is turbulucid have extensive doctrings! Use them to get more information.
In [11]:
?plot_field()
Let us now make a vector plot. This can be done using the plot_vectors() function. We need to pass the Case object and either the name of a vector field present in the case or an array with several column-components. This alternative is implemented for all the plot functions, including plot_field that we looked at above.
In [12]:
a = plot_vectors(case, 'UMean')
This is not very satisfying. Creating a good vector plot usually needs some tweeking of additional keyword arguments. By default the length of the arrows is proportional to the magnitude of the vector, but this can be changed by the normalize keyword.
In [13]:
plot_vectors(case, 'UMean', normalize=True)
Out[13]:
This, perhaps, looks slightly better, but the arrows are obviously too small. This can be controlled by the scale keyword argument. The behaviour of scale is not intuitive: the larger the value of scale the smaller the vectors. This is how the parameter behaves in pyplot.quiver, which is the function that plot_vectors uses to create the plot, so it was decided to leave that as it is. Typically several iterations are needed to get a good value for scale.
Let us make the arrows larger and zoom in on the recirculation zone. Additionally, let us scale the axes.
In [14]:
plot_vectors(case, 'UMean', normalize=True, scaleX=h, scaleY=h, scale=60)
plt.xlim([-1, 6])
plt.ylim([-1, 1])
Out[14]:
OK, this is looking better, we can clearly see two recirculation boubles. The issue is, however, that there are too much arrows. By default, an arrow will be drawn at each cell center, which means the density of arrows is proportional to the density of the mesh. To avoid this, plot_vectors supports another mechanism for positioning the arrows. Instead of using each cell center it is possible to use a Cartesian grid of prescribed resolution to sample the values of the vector field and then plot arrows only at the locations of the nodes of the grid. To do this the sampleByPlane keyword argument should be set to True and the keyword argument planeResolution should provide the amount of nodes in the Cartesian grid in x and y, respectively.
To illustrate the concept let us use a grid with 10 rows and 20 columns.
In [15]:
plot_vectors(case, 'UMean', scaleX=h, scaleY=h, sampleByPlane=True, planeResolution=[10, 20])
Out[15]:
OK, let's now increase the resolution of the mesh, tweak scale, and zoom in on the recirculation zone.
In [16]:
plot_vectors(case, 'UMean', scaleX=h, scaleY=h, normalize=True, sampleByPlane=True, planeResolution=[70, 150], scale=50)
plt.xlim([-1, 6])
plt.ylim([-1, 1])
Out[16]:
Looking good! Let us finally add some color to the arrows. One possibility is to use a single color for all the arrows, this can be done using the color keyword argument, similar to how this is done in e.g. plt.plot() and similar function. It is also possible to color the arrows by the value of field present in a case. In this case, the colorField argument should be used and set to the name of the scalar field that is to be used for coloring. You can also pass an array of values directly, but it is somewhat tricky, when sampling by plane is also used, since the array size and form has to comply to that of the sampling plane. Let us color the vectors using the velocity magnitude field.
In [17]:
plot_vectors(case, 'UMean', colorField='magUMean', scaleX=h, scaleY=h, normalize=True, sampleByPlane=True, planeResolution=[70,150], scale=50)
plt.xlim([-1, 6])
plt.ylim([-1, 1])
Out[17]:
Let us now plot some streamlines. The corresponding function is called plot_streamlines. As with the other plotting functions, we need to pass a Case object and either a field name or an array.
In [18]:
plot_streamlines(case, 'UMean')
Out[18]:
This looks quite nice already. Under the hood the function uses pyplot.streamplot. Unfortunately, pyplot.streamplot does not work with unstructured datasets, and requires data stored on a Cartesian grid. For this reason, plot_streamlines always uses the approach used by plot_vectors when the sampleByPlane keyword argrument is set to True. The plane resolution is set by default to be 50 rows and 50 columns. This, of course, can be adjusted.
As with the vectors, we can add color.
In [19]:
plot_streamlines(case, 'UMean', colorField='magUMean', planeResolution=[300,300], scaleX=h, scaleY=h)
Quite a lot of cutomisation is possible by passing additional keyword arguments, for instance the locations of the seed points for the streamlines. Overall, however, pyplot.streamplot is not the ideal function to rely on and in the future it will probably be attempted to use VTK to compute the streamlines.
Two additionaly plotting functions are probided by the library. The first one, plot_boundary simply plots the boundary of the domain.
In [20]:
plot_boundaries(case, color="Red")
Out[20]:
The second, add_colorbar, adds a colorbar to the plot. This can be useful if more control over the colorbar is needed. By default plot_field adds a colorbar on its own, but this can be avoided by setting the colorbar keyword argument to False.
As it was mentioned in the beginning, the most general way of extracting values from specific locations is working with the cellCentres attribute of the Case class.
However, turbulucid defines several functions that cover commonly needed data extraction procedures.
One of the is profile_along_line which alows to extract a linear profile between two points in the dataset.
In [21]:
point1 = (-0.05, -1)
point2 = (-0.05, 1)
points, data = profile_along_line(case, point1, point2)
Let us explore what has been returned by the function. Firstly, note that we have two returned objects. The first one, assigned to points is a 1D array which holds the locations, along the line, where the sampled data is located.
In [22]:
points
Out[22]:
Note that the first value is 1. This reflects the fact that the line was defined as a vertical segment between y=-1 and y=1, and the first datapoint we have is located at y=0, which then corresponds to a distance of 1 in y from point1.
Quite often it is convenient to subtract the length of the segment, where no data resides, so that the first coordinate is always 0. This can be done with the correctDistance keyword.
In [23]:
points, data = profile_along_line(case, point1, point2, correctDistance=True)
points
Out[23]:
Let us now consider the second returned object. It is a dictionary with the same keys as the Case object, i.e the names of the fields present in the case. Similarly, the values for each key is the values of the corresponding field at locations along the line.
In [24]:
print(data.keys())
data['magUMean']
Out[24]:
Now that we have the data we can use Python's plotting routines to make line-plots.
In [25]:
plt.plot(points, data['pMean'])
Out[25]:
Combining with ploting function we can do more fancy things
In [26]:
plot_boundaries(case, scaleX=h, scaleY=h)
Uref = data['UMean'][-1 ,0]
point1, point2 = (-0.05, 0), (-0.05, 1)
points, data = profile_along_line(case, point1, point2)
plt.plot(data['UMean'][:,0]/Uref - 0.05/h, points/h, 'C1')
point1, point2 = (0, 0), (0, 1)
points, data = profile_along_line(case, point1, point2)
plt.plot(data['UMean'][:,0]/Uref, points/h, 'C1')
point1, point2 = (0.1, -h), (0.1, 1)
points, data = profile_along_line(case, point1, point2)
plt.plot(data['UMean'][:,0]/Uref + 0.1/h, points/h-1, 'C1')
Out[26]:
As it was briefly mentioned in the beginning, turbulucid supports reading data from a multiblock dataset (.vtm), where one block corresponds to the internal field and the other corresponds to defined named boudaries. Functions are provided for accessing both the data on the boundaries, but also from internal cells which are on the boundary, i.e. have one edge that is part of the boundary.
If a .vtk or .vtu file is read, one boundary, called boundary, is defined. Too see what boundaries are present the boundaries attribute of the Case object can be used.
In [27]:
case.boundaries
Out[27]:
Too access data on the boundary, the boundary_data memeber function of the Case object can be used. The name of the boundary is passed as argument.
In [28]:
points, data = case.boundary_data('inletB')
The structure of the returned values is similar to profile_along_line. The first returned object is an array, defining the location of the edge centres of the boundary edges.
In [29]:
points[:10, :]
Out[29]:
The second returned object is the dictionary with the data.
Similarly the boundary_cell_data member function returns the values in the boundary-adjacent cells. The order of the values is the same is that returned by boundary_data.
Another function in turbulucid, dist can be used to compute the distance between the cell centers of the boundary-adjacent cells and a given boundary.
In [30]:
dist(case, 'stepB')
Out[30]:
In this case all the distances are the same due to the geometry being meshed with rectangles.
This covers most of the functionality in turbulucid! More things are available, to explore them look at the code reference!