FunkySetFieldsによるcell zone全体への値の設定方法

openFOAMで使用するメッシュの一部をcell zoneとして陽に定義しておけば、後からcell zone内の変数の値を一括で変更することができる。
これはスカラーの移流拡散やVOFを用いた計算の初期値の設定の際などに応用がきく。

ただ、この文法がいまいちよくわからず、日本語で上記のことを説明している文章が見つからなかったので忘備録として記述する。

ファイル自体に書いておく方法もあるが、これは結構説明してくれている。
http://www.ofwikija.org/images/Swak4Foam%E3%82%92%E4%BD%BF%E3%81%A3%E3%81%A6%E3%81%BF%E3%81%BE%E3%81%97%E3%81%9F%282%29.pdf

コマンドラインからも変更できて便利なので、そっちを書いておく。
なお、swak4foamは必須である。web上の資料もしくは以前のエントリを参照のこと。下記の通り。

funkySetFields -filed 設定したい変数名 -condition "zone(設定したいzone名)" -expression 値 -time 設定したい時間(0なら不要)


平たく言えば、オプションの一つであるconditionに上記のように記述すればよい。
conditionについて、数式を使った設定方法の記述は多いのだが、やりたいことはもっと単純だよねってことで…。

Vsphere-1.7.8 (VCAD)のインストール

Vsphereと言えば、vmwareが出してるソフトウェアの一つだけれども、実は同名のソフトウェアを理研が公開している(要登録)。
VCAD システム研究プログラム

純国産コードのみでのCAEの統合開発環境の構築を試みたVCADシステム、その中核を担うソフトウェアだが、ちょっとググった限りではそれほど広く普及しているとは言い難い。
そんな中で諸事情により本ソフトウェアをローカル環境にインストールする必要に迫られた。

理研の上記リンクからVsphereをダウンロードして展開し、

./configure --prefix=$HOME/Vsphere
make

をまずやってみた。これに失敗。
エラーを見ると、SKLTiming.hの中でfprintfが定義されていないと言われる。いやいや、stdio.hにあるだろ…。
不信感にかられつつSKLTiming.hの中を見るとstdio.hのインクルードがない。別のインクルードファイルで間接的に引用しているとかも見たところなさげ。
わけがわからんなーと思いつつ、SKLTiming.hの一番上に#include を書き足してmakeを行う。

これで最後までmakeができたので

make install

今度も問題なし。まあ、とりあえず動くかな…と、思い、サンプルコードを探してみたら見つからない。どうやらそんな軟弱なものは提供していないらしい。


諸事情から、Vsphereを使った物理計算のコードが手元にあったので、Vsphereのprefixを自分のローカル環境に書き換えてコンパイルしてみた。


しかし再度失敗。どうやらg++やgfortranと相性があまりよろしくないらしい。仕方ないのでIntelコンパイラで再度Vsphereをインストールする。ついでにopenmpiも最新版を設定しておいた(前のエントリ)。

./configure --prefix=$HOME/Vsphere --enable-openmpi=$HOME/openmpi CXX=icpc F90=ifort CXXFLAGS=-O3 F90FLAGS=-O3 LDFLAGS=-L/opt/intel/lib/intel64
make
make install

インストール後、手元のコードをコンパイルすると、無事実行ファイルが作成され、mpiが正常に動作することを確認した。

しかしデフォルトで出力されるファイルはVCADシステムオリジナルのバイナリファイルで中身のフォーマットさえ分からない。計算結果の確認さえできない。
どうやらまた別の可視化ソフトV-isioを使わないといけないらしい。面倒。このへんがいまいち流行ってない原因だろうか…。

openmpi-1.10.1のインストール

openmpiの最新版が二日前にリリースされていたので、ローカル環境にインストールを行った。
tar.gzファイルを落として解凍し、
./configure --prefix=$HOME/openmpi-1.10.1
make
make install
で終了。.bashrcに
export PATH=$HOME/openmpi-1.10.1/bin:$PATH
export MANPATH=$HOME/openmpi-1.10.1/share/man:$MANPATH

と書いて終了。気軽なもんです。includeとlibのパスはmakefileに書くので無視。


2016/04追記
bigvalley.hatenablog.com

C++によるHDF5からの入力

日曜日の全てをかけてHDF5とXDMFの使い方を調べてきたが、これでやっと一段落だろうか。
C++によるHDF5からのデータの入力方法について、公式サイトのsampleがこちら。

HDF5 C++ API: readdata.cpp

いろいろとsampleとしての処理を行っているためか、ぱっと見で把握できなかった。
そこで上を参考にして、ただdouble型の2次元配列を読み込むためだけのコードを作成してみた。
下記に示す。templateを使用すれば任意の型に対応できると思う。
(追記(11/3):boostのmulti_arrayを導入して一部簡略化した)

#include<iostream>
#include<H5Cpp.h>
#include"boost/multi_array.hpp"

using namespace std;

int main(){

	H5::H5File file("SDS.h5", H5F_ACC_RDONLY);
	H5::DataSet dataset = file.openDataSet("doubleArray2"); 
	H5::DataSpace dataspace = dataset.getSpace();

	int rank = dataspace.getSimpleExtentNdims();

	hsize_t dims_out[2];
	int ndims = dataspace.getSimpleExtentDims( dims_out, NULL);
	int NX = (unsigned long)(dims_out[0]);
	int NY = (unsigned long)(dims_out[1]);
	cout << "rank " << rank << ", dimensions " << NX << " x " << NY << endl;

	boost::multi_array<double,2> data(boost::extents[NX][NY]);

	dataset.read(data.data(),H5::PredType::NATIVE_DOUBLE);

	for(int i=0;i<NX;i++){
		for(int j=0;j<NY;j++) cout << data[i][j] << " ";
		cout << endl;
	}
}

先の書き込み用フォーマットと併せてオリジナルのクラスを作成すれば、HDF5におけるIO処理はほぼ全ての処理が自動化できる。
ここまでやって初めて、自分でバイナリのファイルを作成するよりはるかに作業の効率化が望めるだろう。週末をこれでつぶした甲斐はあったかな…。

C++によるHDF5への出力

HDF5のC++での使い方は下記にまとめられている。
HDF5 C++ API: Main Page
かなりわかりにくいけど、一応下記を見れば全てわかるのだろう、わかる人には。今の私には時間がかかりすぎるのでパス。


とにかく他のサンプルコードも見て、わかったこと。
HDF5への出力において一番面倒なことは、出力フォーマットとして一次元のベクトルデータにしか対応していない点である。
数値計算はその多くが行列(テンソル)演算であり、可視性を守るためには速度を犠牲にしつつも多次元配列の形でテンソルを定義したいところである。


この対処法を調べたところ、下記サイトでサンプルコードを公開しているこれまでとは違う神を発見。
今度はOxfordの研究者らしい。
HDF5 C++ interface: writing dynamic 2D arrays - Stack Overflow

C++の拡張ライブラリであるboostに実装されているmultiarrayを使用することで、一次元配列をあたかも多次元配列であるかのように処理しつつ、hdf5への出力では、しれっと一次元配列として入力している、と、見ればいいはず。

実際に試してみる。do_write_hdf5関数の中のvectorにstd::を付け忘れているけど、修正したら問題なく動く。multiarray型の任意の型の多次元配列に対応している。超便利。


vector型は使いやすくて好きだったけど、この分だと早々にboostのmultiarrayへ乗り換えた方がよさげ。
push_backなど、vector型でないとできないような操作はvectorで行うとして、普通の固定長動的配列とみなせる配列は全てmultiarrayで定義したほうがいいかもしれない。

しかしboostの使用はやはりC++userには避けられないのか…。

XDMFによる時系列データの可視化

XDMFファイルでは時系列データの可視化手順について2つのオプションがある。

1.VTK同様、XDMFを連番で出力する(ex. hogehoge_0000.xmf, hogehoge_0001.xmf...)

2 一つのXDMFファイルにまとめる。
全ての時刻歴のデータを一つのXDMFファイルに書き込まないといけなくなるので、連番のファイルが数百等々になると、手軽に可視化できなくなりそう。計算結果の可視化手順が煩雑になるのは何より避けたい。
Xdmf2 Model and Format Archive - XdmfWeb
やり方自体は一応公式サイトにも載っている。
後は前のエントリで紹介した神が、やはりサンプルを書いてくれている。
[Paraview] reading timeseries with XDMFReader?
Duke大のvisulaizationグループの研究者のよう。しかし上のログが2010年のログであることが恐ろしい。最先端から見れば自分は5年以上遅れている。


とにかく、一つの時間でのデータをHDF5により一つのデータにまとめ、その都度XDMFファイルを時間歴に合わせて作成するのが、効率と手軽さを考えると一番楽かな。

Para-view(or Visit)による非構造格子の可視化のためのXDMFファイル(.xmf)の作成

任意の非構造格子の可視化ファイルを下記に示す。
[visit-users] using hdf5 files with visit (XDMF help)
ここで提供されているサンプルのうち一つに変数を付け加えただけである。

可視化に最低限必要な情報は節点数、要素数、節点座標、各要素が持つ節点の番号のみであり、vtkやplyファイルと同様である。
ただ、可視化できるデータの量は格段に増えている。一度のファイルでほぼ無制限に節点や要素が持つ変数を設置できる。これは普通のvtkファイルだとできなかった(はず)。
変数もスカラー、ベクトルはおろか、テンソルにも対応している。恐ろしい子である。



下記はvtkなどの可視化フォーマットに慣れている人ならばすぐに理解できると思われる。
そうでなくても、xml形式で各変数の意味を説明しているので、下記をひな型に自分の用途に合わせて編集するのは簡単だろう。
Dimensionを定義しているため、数値の改行は本来必要ない。

<?xml version="1.0" ?>
<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>

<Xdmf>
    <Domain>
        <Grid Name="Mixed">
            <Topology Type="Mixed" NumberOfElements="3" >
                <DataItem Format="XML" DataType="Int" Dimensions="19">
		6   0 1 2 6
		5   3 4 5 6
                9   7 8 9 10 11 12 13 14
                </DataItem>
            </Topology>
            <Geometry Type="XYZ">
                <DataItem Format="XML" Dimensions="15 3">
                0.0    0.0    0.0
                1.0    0.0    0.0
                1.0    1.0    0.0
                0.0    0.0    2.0
                1.0    0.0    2.0
                1.0    1.0    2.0
                0.0    1.0    2.0
                0.0    0.0    3.0
                1.0    0.0    3.0
                1.0    1.0    3.0
                0.0    1.0    3.0
                0.0    0.0    4.0
                1.0    0.0    4.0
                1.0    1.0    4.0
                0.0    1.0    4.0
                </DataItem>
            </Geometry>
            <Attribute Name="NodeScalar" Center="Node">
                <DataItem Format="XML" Dimensions="15">
		 100 200 300 400 500
		 600 600 700 800 900
	     1000 1100 1200 1300 1400
                </DataItem>
	</Attribute>
        <Attribute Name="CellScalar" Center="Cell">
            <DataItem Format="XML" Dimensions="3">
	    	1 2 3 
            </DataItem>
        </Attribute>
        <Attribute Name="NodeVector" AttributeType="Vector" Center="Node">
            <DataItem Format="XML" Dimensions="15 3">
	    	1.0 0.0 0.0 
	    	2.0 0.0 0.0  
	    	3.0 0.0 0.0 
	    	4.0 0.0 0.0
	    	5.0 0.0 0.0
	    	6.0 0.0 0.0
	    	6.0 0.0 0.0
	        7.0 0.0 0.0 
	    	8.0 0.0 0.0
	    	9.0 0.0 0.0
	    	10.0 0.0 0.0
	        11.0 0.0 0.0
	    	12.0 0.0 0.0
	    	13.0 0.0 0.0
	    	14.0 0.0 0.0
            </DataItem>
	</Attribute>
        <Attribute Name="CellVector" AttributeType="Vector" Center="Cell">
            <DataItem Format="XML" Dimensions="3 3">
	    	1 0 0
	    	2 0 0
	    	3 0 0
            </DataItem>
        </Attribute>
        <Attribute Name="NodeTensor" AttributeType="Tensor" Center="Node">
            <DataItem Format="XML" Dimensions="15 9">
	    	1.0 2.0 3.0 1.0 2.0 3.0 1.0 2.0 3.0
	    	2.0 1.0 3.0 1.0 3.0 2.0 1.0 3.0 2.0
	    	3.0 4.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0
	    	4.0 0.0 0.0 4.0 5.0 4.0 4.0 5.0 4.0
	    	5.0 0.0 0.0 1.0 2.0 3.0 1.0 2.0 3.0
	    	6.0 0.0 0.0 1.0 3.0 2.0 1.0 3.0 2.0
	    	6.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
	        7.0 0.0 0.0 4.0 5.0 4.0 4.0 5.0 4.0
	    	8.0 0.0 0.0 1.0 2.0 3.0 1.0 2.0 3.0
	    	9.0 0.0 0.0 1.0 3.0 2.0 1.0 3.0 2.0
	    	1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
	        1.0 0.0 0.0 4.0 5.0 4.0 4.0 5.0 4.0
	    	1.0 0.0 0.0 1.0 2.0 3.0 1.0 2.0 3.0
	    	1.0 0.0 0.0 1.0 3.0 2.0 1.0 3.0 2.0
	    	1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            </DataItem>      
	</Attribute>
        <Attribute Name="CellTensor" AttributeType="Tensor" Center="Cell">
            <DataItem Format="XML" Dimensions="3 9">
	    	1 0 0 1.0 2.0 3.0 1.0 2.0 3.0
	    	2 0 0 1.0 3.0 2.0 1.0 3.0 2.0
	    	3 0 0 0.0 0.0 0.0 0.0 0.0 0.0
            </DataItem>
       </Attribute>
        <Attribute Name="NodeSymTensor" AttributeType="Tensor6" Center="Node">
            <DataItem Format="XML" Dimensions="15 9">
	    	1.0 2.0 3.0 1.0 2.0 3.0
	    	2.0 1.0 3.0 1.0 3.0 2.0
	    	3.0 4.0 1.0 0.0 0.0 0.0
	    	4.0 0.0 0.0 4.0 5.0 4.0
	    	5.0 0.0 0.0 1.0 2.0 3.0
	    	6.0 0.0 0.0 1.0 3.0 2.0
	    	6.0 0.0 0.0 0.0 0.0 0.0
	        7.0 0.0 0.0 4.0 5.0 4.0
	    	8.0 0.0 0.0 1.0 2.0 3.0
	    	9.0 0.0 0.0 1.0 3.0 2.0
	    	1.0 0.0 0.0 0.0 0.0 0.0
	        1.0 0.0 0.0 4.0 5.0 4.0
	    	1.0 0.0 0.0 1.0 2.0 3.0
	    	1.0 0.0 0.0 1.0 3.0 2.0
	    	1.0 0.0 0.0 0.0 0.0 0.0
            </DataItem>      
	</Attribute>
        <Attribute Name="CellSymTensor" AttributeType="Tensor6" Center="Cell">
            <DataItem Format="XML" Dimensions="3 9">
	    	1 0 0 1.0 2.0 3.0
	    	2 0 0 1.0 3.0 2.0
	    	3 0 0 0.0 0.0 0.0
            </DataItem>
       </Attribute>
        </Grid>
    </Domain>
</Xdmf>

上のxmlを可視化したものが以下。二階のテンソルを入力しているので、次元が9つある。

f:id:bigvalley:20151101175758p:plain

要素の定義として、要素の種類を番号で設定する必要がある(6:四面体,5:四角形,9:六面体)が、下記の公式サイトにはその資料が見つからなかった。
XdmfWeb

調べたところ、XDMFをさらに拡張したフォーマットを提供しているグループが説明書きを書いてくれていた。
Separated Variables Representation Visualisation | Reduced Order Modeling


上記サイトで提供しているpdfファイルに細かく書いてある。

なお、データ数が増加する場合はxml形式に変数をべた書きするのではなく、
hdf5ファイルなどの外部ファイルにバイナリ形式でまとめたほうが、データの扱いとしてスマートだし可視化速度も速いと思う。

上記は簡単な例だけども、xmlの書き方次第でマルチドメインにも簡単に対応できる。