OpenFoam is a popular open source code for computational fluid dynamics (CFD). Although it contains various helpful postpocessing modules in command line such as postProcess, it is still designed for convenience but not flexibility. For example, it only provides operations that are very common in the context of fluid mechanics or vector mathematics and it hides the details of operation such as the numerical scheme to approximate the derivatives. Most of the time, OpenFoam saves the data in folder named by the current time and in each folder contains a special OpenFoam format-txt like data, which is also designed for convenience such that one can directly read the result in the field data. However, if one wants to manipulate the data in a more flexible sense in the modern data-driven era, a Python-script-driven manipulation of data is extremely favorable. Also, to avoid dealing with the mesh in the script, if would be great if one can simply add the modified field on the original mesh.
Fortunately, with the help of an awesome Python package on Github: currently Pyvista previously vtkInterface , originally by Alex Kaszynski, one can easily leverage the powerful libraries in Python environment to postprocessing traditional, mature, standard and specialized scientific computing data and immediately put them back in to again, leverage the existing powerful visualization software in scientific computing community.
Using pip to install vtkInterface
Note that it supports well in Python 3.5+.
sudo pip install pyvista
Tutorial: 2D Flow past cylinder
untarthe .tar file
tar -zxvf vortex_shedding.tar.gz ./
c1directory for running a standard case
cd c1 blockMesh checkMesh icoFoam > log &
Convert OpenFoam default format to VTK
Using Python to manipulate VTK data
import pyvista as vtki import numpy as np ## grid is the central object in VTK where every field is added on to grid grid = vtki.UnstructuredGrid('./VTK/c1_1000.vtk') ## point-wise information of geometry is contained print grid.points ## get a dictionary contains all cell/point information print grid.cell_arrays # note that cell-based and point-based are in different size print grid.point_arrays # ## get a field in numpy array p_cell = grid.cell_arrays['p'] ## create a new cell field of pressure^2 p2_cell = p_cell**2 grid._add_cell_scalar(p2_cell, 'p2') ## remember to save the modified vtk grid.save('./VTK/c1_1000_shaowu.vtk')
Visualize the new field in ParaView