OpenFOAMのpostProcessコマンドでVTK出力

概要

OpenFOAM 4以前ではsampleコマンドで分布やコンターのデータを取得できてましたが，ver 4以降ではsampleコマンドに代わってpostProcessで処理するようになりました。postProcessコマンドでのデータ出力の手法についてまとめます。

ここではコンターをかくことを想定して，vtkデータを出力します。LES等の解析でデータ量が大きいと，paraFoamで開くのに時間が非常にかかります。またさらにそこからコンターを書こうとしても時間がかかります。そこでvtkではじめからコンター用のデータを出力しておいて，可視化にかかる時間節約を行います。

手順

foamGet

公式の説明（ https://cfd.direct/openfoam/user-guide/v6-graphs-monitoring/ ）にあるとおりですが，はじめにfoamGetコマンドでサンプルを作るのが手っ取り早いです。コンター出力にはオプションでsurfacesを指定します。

$ foamGet surfaces

これでsystem以下にsurfacesというファイルが作成されます。
caseDicts/postProcessing/visualization/surfaces.cfg（特に設定を変えていない場合はOpenFOAM内のetc以下にcaseDictsがあります）がリファレンスとして使われます。なので$cuttingPlaneなどはsurfaces.cfgですでに定義されている値に置き換わります。

cuttingPlaneがサンプルファイルでのtypeになっています。type planeの方がデータを節約できるので，以下のように修正します。

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  6
     \\/     M anipulation  |
-------------------------------------------------------------------------------
Description
    Writes out surface files with interpolated field data in VTK format, e.g.
    cutting planes, iso-surfaces and patch boundary surfaces.

    This file includes a selection of example surfaces, each of which the user
    should configure and/or remove.

\*---------------------------------------------------------------------------*/

# includeEtc "caseDicts/postProcessing/visualization/surfaces.cfg"

fields       (T U p);

surfaces
(
    yz
    {
        type                 plane;
        planeType            pointAndNormal;    // always triangulated
        pointAndNormalDict
        {
            basePoint       (0 0 0);
            normalVector    (1 0 0);
        }
    }
);

// ************************************************************************* //

• fieldsが出力する変数
• OpenFOAM 4より前だと，type planeを指定するときににplaneTypeのオプションは必要なかったみたいですが，今は必要です。

postProcess

surfacesのファイルが修正できたら，postProcessを実行します。オプションにsurfacesを与えます。

$ postProcess -func surfaces

うまく行けば，postProcessing/surfaces以下にデータが出力されます。

caseDicts/postProcessing/visualization/surfaces.cfgで出力形式がvtkになっているので，vtkファイルが作成されます。gnuplotでプロットするならsurfaceFormat raw;にします。（多分。確かめてないけど）ただし，pm3d使ってコンターを書くときにgnuplotだと，データに適当なところで空白行がないとpm3dの処理で失敗します。データサイズが大きいと，手で修正するのもつらいので，gnuplotによるコンター作図はやめました。vtk出力にして，matplotlibでの作図に乗り換えました。

蛇足：OpenFOAMで出力されるVTKファイルの精度

これでVTKファイルが出力できました。ここで若干問題が発生するときがあります。というか問題が出ました。
多成分気体の計算をしてますが，この解析では各ガス種の質量分率が変数として出てくるので，そのデータをVTKを作りました。VTKデータ自体は特に問題なく生成されるのに，いざparaviewでVTKを見てみると，値が1e33とかになっているところがあります。分率なので0から1の値にしかならないはずなのにとてもおかしい。

結論からいうと，VTKでのcell dataの精度がfloatなのが原因でした。0に非常に近い値として1e-39オーダーのデータが含まれていたため，単精度浮動小数のレンジからはみ出すので，paraviewで可視化するとおかしなことになるという結果でした。cell dataの精度をdoubleにしたら無事解決。OpenFOAMの計算では倍精度にしてますが，可視化用のデータ出力の精度は元の計算結果とリンクしてないっぽいです。

単精度浮動小数点で表現できる最小の正の数が1e-38オーダーだってのを初めて知りました。https://ja.wikipedia.org/wiki/単精度浮動小数点数