[PYTHON] Try drawing contour plots using matplotlib in OpenFOAM
Introduction
This article is about tips of OpenFOAM, which is open source software for fluid analysis. In OpenFOAM, the calculation result is drawn using ParaView as post-processing (visualization), but basically the drawn image is in raster format and is two-dimensional. Contour lines (contour) in the cross section cannot be output in vector format. Also, when the cross section and variables to be visualized are decided, the work of starting ParaView and drawing becomes troublesome. Therefore, this article summarizes how to semi-automatically draw a 2D section contour and output it in vector format using matplotlib, a visualization library for python.
Environment
The software required this time is OpenFOAM and Python. The version of OpenFOAM has recently been major updated to 4.0 or v1606, but this article focuses on 2.4.x. For Python, 2.7 is the target. Numpy, matplotlib, vtk (python-vtk), tkinter (python-tk) are required as Python packages, so if not, please install them additionally.
After building the environment, check the operation with the terminal.
OpenFOAM
$ foamInstallationTest
$ mkdir -p $FOAM_RUN
$ run
$ cp -r $FOAM_TUTORIALS/incompressible/icoFoam/cavity ./
$ cd cavity
$ blockMesh
$ icoFoam
foamInstallationTest checks the basic environment settings of the installed OpenFOAM. In the following commands, the cavity flow, which is the basic tutorial of OpenFOAM, is executed, and if there are no errors up to this point, it is OK.
Python
$ python
>>> import numpy
>>> import matplotlib
>>> import vtk
>>> import Tkinter
>>> quit()
If each package (module) is installed, nothing will be displayed and the prompt will be returned normally. ʻImportError: No module If you get the named xxxx` error, reinsert the appropriate package.
Acquisition of 2D section data (vtk)
First, obtain the 2D cross-section data you want to draw from the OpenFOAM analysis results in vtk format. You can read OpenFOAM data with ParaView and output cross-section data with Slice, but the purpose of this time is to automatically draw contour maps with commands or scripts, so ParaView that requires GUI operation is Not used. Here, we will use the OpenFOAM utility sample (which is integrated into postProcess in OpenFOAM 4.0 and later).
Create cross-section data using the results of the cavity tutorial performed in the environment setting check. First, in order to use sample, copy and edit the template of sampleDict, which is the configuration file of this command.
$ run
$ cd cavity
$ cp $FOAM_UTILITIES/postProcessing/sampling/sample/sampleDict ./system
$
Various sampling methods including usage comments are written in the original template file, but this time only the $ x sampleDict is simple as follows. It becomes.
setFormat raw;
surfaceFormat vtk;
formatOptions
{
}
interpolationScheme cellPoint;
fields
(
p
U
);
sets
(
);
surfaces
(
interpolatedPlane
{
type plane;
coordinateSystem
{
origin (0.05 0.05 0.005);
coordinateRotation
{
type axesRotation;
e1 (1 0 0);
e2 (0 1 0);
}
}
basePoint (0 0 0);
normalVector (0 0 1);
interpolate true;
}
);
Since sets is used when acquiring one-dimensional data such as straight lines, leave it blank this time. There are several ways to specify surfaces, but this time we will use the one obtained by interpolation.
After editing sampleDict, execute the sample command to generate 2D section data. Only the data of the last time (0.5) is sampled with the option -latestTime.
$ sample -latestTime
When the command is executed, a directory called postProcessing / surfaces will be created in the current directory, and vtk data at the specified time will be output.
Python script preparation
The format of the vtk file created by the above sample utility is the so-called legacy version 2.0. In Python, if you can read vtk format data and pass the data to the visualization library matplotlib, you can draw freely. You can write your own script, but Alexey Matveichev wrote [Plotting VTK files with matplotlib and python-vtk](http://matveichev.blogspot.jp/2014/03/plotting-vtk-files-with-matplotlib] -and.html) is written, so I will use the published script. The script published in the article has a vtk version of python-vtk of 5.x series, so if you build an environment with the new 6.x series, it is linked in the article [Script published on github](https:: Please use //gist.github.com/mrklein/44a392a01fa3e0181972). Here, we will proceed assuming that the 5.x system is used.
Copy the original script and place it in the OpenFOAM parsing directory (let's name it plotUx.py) and modify it a bit.
plotUx.py
#!/usr/bin/env python
import sys # new
import os # new
import numpy as N
import vtk
import matplotlib.pyplot as P
def load_velocity(filename):
if not os.path.exists(filename): return None
reader = vtk.vtkPolyDataReader()
reader.SetFileName(filename)
reader.ReadAllVectorsOn()
reader.Update()
data = reader.GetOutput()
cells = data.GetPolys()
triangles = cells.GetData()
points = data.GetPoints()
point_data = data.GetPointData()
Udata = point_data.GetArray(0)
ntri = triangles.GetNumberOfTuples()/4
npts = points.GetNumberOfPoints()
nvls = Udata.GetNumberOfTuples()
tri = N.zeros((ntri, 3))
x = N.zeros(npts)
y = N.zeros(npts)
ux = N.zeros(nvls)
uy = N.zeros(nvls)
for i in xrange(0, ntri):
tri[i, 0] = triangles.GetTuple(4*i + 1)[0]
tri[i, 1] = triangles.GetTuple(4*i + 2)[0]
tri[i, 2] = triangles.GetTuple(4*i + 3)[0]
for i in xrange(npts): # modified
pt = points.GetPoint(i)
x[i] = pt[0]
y[i] = pt[1]
for i in xrange(0, nvls): # modified
U = Udata.GetTuple(i)
ux[i] = U[0]
uy[i] = U[1]
return (x, y, tri, ux, uy)
P.clf()
fn = sys.argv[1] # new
levels = N.linspace(-1.0, 1.0, 21) # modified
x, y, tri, ux, uy = load_velocity(fn) # modified
P.tricontour(x, y, tri, ux, levels, linestyles='-',
colors='black', linewidths=0.5)
P.tricontourf(x, y, tri, ux, levels) # modified
P.xlim([0, 0.1]) # new
P.ylim([0, 0.1]) # new
P.gca().set_aspect('equal')
P.minorticks_on()
P.gca().set_xticklabels([])
P.gca().set_yticklabels([])
P.title('$U_x$') # modified
P.savefig("Ux.png ", dpi=300, bbox_inches='tight')
P.savefig("Ux.pdf", bbox_inches='tight')
P.show() # new
The major changes from the original script are that the vtk file to be read is given as an argument using sys.argv, and finally the figure is displayed on the screen withshow (). It is a point.
Execution of Python script (contour diagram of x component of velocity)
The prepared script is for reading vector data, so let's execute it on the sample speed data.
$ python plotUx.py postProcessing/surfaces/0.5/U_interpolatedPlane.vtk
If the contour map is displayed in another window without any error, it is successful. To exit, close the displayed window. Since it is a script that also outputs png and pdf, you can confirm that those image files are created.
Vorticity contour map
So far, I have drawn the $ x $ component of velocity, but I will also draw other physical quantities.
Vorticity calculation
At the stage when the cavity calculation of OpenFOAM is completed, velocity and pressure are output as the main physical quantities. In fluid analysis, the vorticity is used to observe the vortices existing in the analysis area, so let's also visualize the vorticity. First, run the OpenFOAM standard vorticity calculation utility vorticity.
$ vorticity -latestTime
Since the option -latestTime is added, the vorticity is calculated only for the last time. The vorticity data is output under the 0.5 directory.
Cross-section data sampling
Edit the sampleDict prepared for speed and add the field variables to be sampled as follows.
...
fields
(
p
U
vorticity // new
);
...
Execute the sample command in the same way as for speed.
$ sample -latestTime
After execution, vtk data of vorticity (vorticity) will be output under the directory postProcessing / surfaces / 0.5.
Preparation of Python script for vorticity
For a two-dimensional cavity flow, the vorticity has a $ z $ directional component perpendicular to the $ x plotWz.py.
plotWz.py
#!/usr/bin/env python
import sys
import os
import numpy as N
import vtk
import matplotlib.pyplot as P
def load_vorticity(filename): # modified
if not os.path.exists(filename): return None
reader = vtk.vtkPolyDataReader()
reader.SetFileName(filename)
reader.ReadAllVectorsOn()
reader.Update()
data = reader.GetOutput()
cells = data.GetPolys()
triangles = cells.GetData()
points = data.GetPoints()
point_data = data.GetPointData()
Wdata = point_data.GetArray(0) # modified
ntri = triangles.GetNumberOfTuples()/4
npts = points.GetNumberOfPoints()
nvls = Wdata.GetNumberOfTuples() # modified
tri = N.zeros((ntri, 3))
x = N.zeros(npts)
y = N.zeros(npts)
wx = N.zeros(nvls) # modified
wy = N.zeros(nvls) # modified
wz = N.zeros(nvls) # new
for i in xrange(0, ntri):
tri[i, 0] = triangles.GetTuple(4*i + 1)[0]
tri[i, 1] = triangles.GetTuple(4*i + 2)[0]
tri[i, 2] = triangles.GetTuple(4*i + 3)[0]
for i in xrange(npts):
pt = points.GetPoint(i)
x[i] = pt[0]
y[i] = pt[1]
for i in xrange(0, nvls):
W = Wdata.GetTuple(i) # modified
wx[i] = W[0] # modified
wy[i] = W[1] # modified
wz[i] = W[2] # new
return (x, y, tri, wx, wy, wz) # new
P.clf()
fn = sys.argv[1]
levels = N.linspace(-260, 260, 16) # modified
x, y, tri, wx, wy, wz = load_vorticity(fn) # modified
P.tricontour(x, y, tri, wz, levels, linestyles='-',
colors='black', linewidths=0.5) # modified
P.tricontourf(x, y, tri, wz, levels) # modified
P.xlim([0, 0.1])
P.ylim([0, 0.1])
P.gca().set_aspect('equal')
P.minorticks_on()
P.gca().set_xticklabels([])
P.gca().set_yticklabels([])
P.title('$W_z$') # modified
P.savefig("Wz.png ", dpi=300, bbox_inches='tight') # modified
P.savefig("Wz.pdf", bbox_inches='tight') # modified
P.show()
Execution of Python script (contour diagram of z component of vorticity)
Run the script as you would for speed.
$ python plotWz.py postProcessing/surfaces/0.5/vorticity_interpolatedPlane.vtk
If the contour map is displayed in another window without any error, it is successful.
Contour diagram of pressure
I have drawn the $ x $ component of velocity and the $ z $ component of vorticity, but there is also pressure as the output physical quantity. Since the pressure is a scalar value, you only need to change the value read in the first speed script from two to one. Since the modification of the script is simple, it is omitted here. If you are interested, please create it yourself. Since the pressure (p) is already specified in the sampleDict of OpenFOAM, you can draw a contour map of the pressure by executing the modified script.
Summary
In this article, I introduced how to visualize the calculation result of OpenFOAM as a contour map with a two-dimensional cross section using matplotlib. Depending on the OS you are using, it may be troublesome to build a Python environment, but once the construction is complete, you can change the calculation conditions and automatically generate a visualization diagram. This time, the contour line spacing (levels) is entered manually, but it is also possible to find the maximum / minimum and automatically determine it on Python. In addition, matplotlib can draw various visualization diagrams, so please customize it by yourself along with the functions of Python.
Recommended Posts