OpenFOAMでmatplotlibを用いて等高線図を描いてみる

はじめに

この記事は流体解析のオープンソースソフトウェアであるOpenFOAMのtipsに関する記事です．OpenFOAMでは後処理(可視化)としてParaViewを用いて計算結果の描画を行いますが，基本的に描画された画像はラスター形式であり，二次元断面での等高線(コンター)をベクター形式で出力できません．また，可視化する断面・変数が決まっているときはParaViewを起動して描画する作業が面倒になります．そのため，この記事ではpythonの可視化ライブラリであるmatplotlibを用いて半自動的に二次元断面コンターを描画してベクター形式で出力する方法をまとめています．

環境構築

今回必要とするソフトウェアはOpenFOAMとPythonです．OpenFOAMのバージョンが最近メジャーアップデートされて4.0もしくはv1606になっていますが，この記事では2.4.xを対象とします．またPythonに関しては2.7を対象とします．Pythonパッケージとしてnumpy，matplotlib，vtk (python-vtk)，tkinter (python-tk) を必要としますので無ければ追加インストールしてください．

環境構築が済んだら，端末ターミナルによって動作を確認してみます．

OpenFOAM

$ foamInstallationTest
$ mkdir -p $FOAM_RUN
$ run
$ cp -r $FOAM_TUTORIALS/incompressible/icoFoam/cavity ./
$ cd cavity
$ blockMesh
$ icoFoam

foamInstallationTestはインストールされたOpenFOAMの基本的な環境設定をチェックしています．以降のコマンドではOpenFOAMの基本チュートリアルであるキャビティ流れを実行しており，ここまででエラーが出ていなければOKです．

Python

$ python
>>> import numpy
>>> import matplotlib
>>> import vtk
>>> import Tkinter
>>> quit()

各パッケージ(モジュール)がインストールされていれば何も表示されずに正常にプロンプトが返ってきます．ImportError: No module named xxxxのエラーが返ってきたら，もう一度該当するパッケージを入れ直してください．

二次元断面データ(vtk)の取得

まずOpenFOAMの解析結果から，描画したい二次元断面データをvtk形式で取得します．OpenFOAMのデータをParaViewで読み込んでSliceで断面データを出力してもよいのですが，今回の趣旨はなるべくコマンドorスクリプトで自動的に等高線図を描くことなので，GUI操作が必要なParaViewは用いません．ここではOpenFOAMのユーティリティであるsample(OpenFOAM 4.0以降ではpostProcessに統合されている)を利用します．

環境設定確認のところで実行したキャビティチュートリアルの結果を使って断面データを作成します．まずsampleを使うためにこのコマンドの設定ファイルであるsampleDictの雛形をコピーして編集します．

$ run
$ cd cavity
$ cp $FOAM_UTILITIES/postProcessing/sampling/sample/sampleDict ./system
$ 

元の雛形ファイルには使い方コメントを含めて様々なサンプリング方法が書いてありますが，今回は二次元キャビティの$x$-$y$平面だけをサンプルするので，sampleDictを以下のように単純化します．

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;
    }
);

setsは直線などの一次元データを取得するときに使うものなので今回は空白にします．surfacesの指定の仕方は何通りかありますが，今回は補間によって求めるものを使用します．

sampleDictの編集が終わったら，sampleコマンドを実行して二次元断面データを生成します．オプション-latestTimeを付けて最終時刻(0.5)のデータのみサンプルします．

$ sample -latestTime

コマンドを実行すると，カレントディレクトリにpostProcessing/surfacesというディクレトリができ，指定した時刻でのvtkデータが出力されます．

Pythonスクリプトの準備

上述のsampleユーティリティで作成されたvtkファイルのフォーマットはいわゆるレガシー形式のバージョン2.0です．Pythonではvtk形式のデータを読み込み，可視化ライブラリのmatplotlibにデータを受け渡すことができれば自由に描画できます．自分でスクリプトを書いても良いですが，Alexey MatveichevさんがPlotting VTK files with matplotlib and python-vtkという記事を書かれているので，公開されているスクリプトを利用します．記事で公開されているスクリプトはpython-vtkのvtkバージョンが5.x系のものなので，新しい6.x系で環境構築した場合は記事中でリンクされているgithub公開のスクリプトを利用してください．ここでは5.x系を利用するとして進めます．

元のスクリプトをコピーしてOpenFOAMの解析ディレクトリに置き(ファイル名はplotUx.pyとしておく)，少し修正します．

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

オリジナルのスクリプトからの大きな変更点は，読み込むvtkファイルをsys.argvを用いて引数で与えられるようにした点と，最後にshow()で画面に図が表示されるようにした点です．

Pythonスクリプトの実行 (速度のx成分の等高線図)

準備したスクリプトはベクトルデータ読み込み用なので，サンプルした速度データに対して実行してみます．

$ python plotUx.py postProcessing/surfaces/0.5/U_interpolatedPlane.vtk

エラーなしで別のウィンドウで等高線図が表示されれば成功です．終了するには表示されたウィンドウを閉じます．pngとpdfでも出力するスクリプトになっているので，それらの画像ファイルができていることが確認できます．

渦度の等高線図

ここまでで速度の$x$成分を描画してみましたが，他の物理量も描画してみます．

渦度の計算

OpenFOAMのキャビティ計算を終えた段階では，主な物理量としては速度と圧力が出力されています．流体解析では渦度が解析領域に存在する渦を観察するために用いられるので，渦度も可視化してみます．まずはOpenFOAM標準の渦度計算ユーティリティvorticityを実行します．

$ vorticity -latestTime

オプション-latestTimeを付けたので最終時刻に対してのみ渦度が計算されます．0.5ディレクトリの下にvorticityデータが出力されます．

断面データのサンプリング

速度のときに用意したsampleDictを編集し，サンプルするフィールド変数を以下のように追記します．

...
fields
(
    p
    U
    vorticity  // new
);
...

速度のときと同様にsampleコマンドを実行します．

$ sample -latestTime

実行後に，ディクレトリpostProcessing/surfaces/0.5の下に渦度(vorticity)のvtkデータが出力されます．

渦度用Pythonスクリプトの準備

二次元キャビティ流れの場合，渦度は$x$-$y$平面に垂直な$z$方向成分を持つので，渦度ベクトルの$z$成分を描画する必要があります．そのため，速度$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()

Pythonスクリプトの実行 (渦度のz成分の等高線図)

速度のときと同様にスクリプトを実行します．

$ python plotWz.py postProcessing/surfaces/0.5/vorticity_interpolatedPlane.vtk

エラーなしで別のウィンドウで等高線図が表示されれば成功です．

圧力の等高線図

速度の$x$成分，渦度の$z$成分を描画して来ましたが，出力されている物理量としては圧力もあります．圧力はスカラー値なので，はじめの速度用スクリプトにおいて読み込む値を2つから1つに変更するだけでOKです．スクリプトの修正は単純なのでここでは省略します．興味のある方は自分で作成してみてください．OpenFOAMのsampleDictではすでに圧力(p)を指定しているので，修正したスクリプトを実行すれば圧力の等高線図が描けると思います．

まとめ

この記事ではmatplotlibを利用してOpenFOAMの計算結果を二次元断面で等高線図として可視化する方法を紹介しました．利用しているOSによってはPythonの環境構築が面倒な場合がありますが，構築が済んでしまえば計算条件を変えて可視化図を自動的に生成することが可能です．今回，等高線の間隔(levels)はマニュアルで入力していますが，最大・最小を求めたりしてPython上で自動に決定することも可能です．またmatplotlibは多様な可視化図を描画できますので，Pythonの機能と併せて自分でカスタマイズしてみてください.