Siloとは
LLNLが開発している科学計算用のI/Oライブラリで、
様々なデータ形式用の出力に対応しています。
https://wci.llnl.gov/simulation/computer-codes/silo
同様にLLNLが開発している可視化ソフトウェアViSItがSilo形式を完全サポートしているのが特徴で、
データ出力から可視化までスムーズに行えます。
日本語情報が見つからなかったので、ビルドからC++での利用まで一通りまとめた内容を書いておきます。
ビルド
$ wget https://wci.llnl.gov/content/assets/docs/simulation/computer-codes/silo/silo-4.10.2/silo-4.10.2.tar.gz
$ tar zxvf silo-4.10.2.tar.gz
$ cd silo-4.10.2
$ ./configure --prefix=$HOME/local
$ make
$ make install
prefixはインストールしたい場所に応じて変更する。
C/C++からの利用法
silo.hをインクルードすればOK。
(上記手順でインストールしたなら$HOME/local/include/silo.hにある)
ドキュメント https://wci.llnl.gov/content/assets/docs/simulation/computer-codes/silo/LLNL-SM-654357.pdf がしっかり書いてあるので、利用に困ることなさそう。
また、https://wci.llnl.gov/simulation/computer-codes/silo/examples に特定の形式で出力したい場合のサンプルコードが置いてある。
とりあえず2章に示してあるサンプル通りに書いて出力してみる。
#include <silo.h>
int main(int argc, char* argv[])
{
// Create the Silo file
DBfile* file = DBCreate("sample.silo", DB_CLOBBER, DB_LOCAL, NULL, DB_PDB);
// names of the coordinates
char* coordnames[2] = {"x", "y"};
// dimensions
int dimensions[2];
dimensions[0] = 4;
dimensions[1] = 4;
// the coordinate arrays
float node_x[4] = {-1.1, -0.1, 1.3, 1.7};
float node_y[4] = {-2.4, -1.2, 0.4, 0.8};
// the array of coordinate arrays
float* coordinates[2];
coordinates[0] = node_x;
coordinates[1] = node_y;
DBPutQuadmesh(file, "mesh1", coordnames, coordinates, dimensions, 2, DB_FLOAT, DB_COLLINEAR, NULL);
DBClose(file);
return 0;
}
得られた結果をVisItで可視化した結果がこちら。
sourceにsample.siloを選んで、Add->Mesh->mesh1してactive windowにDrawする。
meshが描画された。便利！！！
三次元の場合
DBPutQuadmeshの定義は
int DBPutQuadmesh (DBfile *dbfile, char const *name,
char const * const coordnames[], void const * const coords[],
int dims[], int ndims, int datatype, int coordtype,
DBoptlist const *optlist);
となっているので、
ndims = 3として、各配列を三次元にすればOK。
#include <silo.h>
int main(int argc, char* argv[])
{
const int dim = 3;
const int nx = 5;
const int ny = 5;
const int nz = 5;
// Create the Silo file
DBfile* file = DBCreate("sample.silo", DB_CLOBBER, DB_LOCAL, NULL, DB_PDB);
// names of the coordinates
char* coordnames[dim] = {"x", "y", "z"};
// dimensions
int dimensions[dim];
dimensions[0] = nx;
dimensions[1] = ny;
dimensions[2] = nz;
// the coordinate arrays
float node_x[nx] = {-2.0, -1.0, 0.0, 1.0, 2.0};
float node_y[ny] = {-2.0, -1.0, 0.0, 1.0, 2.0};
float node_z[nz] = {-2.0, -1.0, 0.0, 1.0, 2.0};
// the array of coordinate arrays
float* coordinates[dim];
coordinates[0] = node_x;
coordinates[1] = node_y;
coordinates[2] = node_z;
DBPutQuadmesh(file, "mesh1", coordnames, coordinates, dimensions, dim, DB_FLOAT, DB_COLLINEAR, NULL);
// Close the Silo file
DBClose(file);
return 0;
}
VisItで出力した図
三次元データの出力
メッシュは出力できたので、対応するデータを出力します。
Quadmeshに対応するデータを入力するのはDBPutQuadvar()で、
その定義は
int DBPutQuadvar (DBfile *dbfile, char const *name,
char const *meshname, int nvars,
char const * const varnames[], void const * const vars[],
int dims[], int ndims, void const * const mixvars[],
int mixlen, int datatype, int centering,
DBoptlist const *optlist);
となっています。
ちょっとハマったのが、vars[]として渡す配列の内容で、
データが1次元の場合でもvars = data;ではなくて、vars[0] = data;という風になってなくてはいけない。
vars:
Array of lengthnvarscontaining pointers to arrays defining the values associated with each subvariable. For true edge- or face-centering (as opposed toDB_EDGECENTcentering whenndimsis 1 andDB_FACECENTcentering whenndimsis 2), each pointer here should point to an array that holdsndimssub-arrays, one for each of the i-, j-, k-oriented edges or i-, j-, k-intercepting faces, respectively. Read the description for more details.
サンプル
#include <silo.h>
int main(int argc, char* argv[])
{
const int dim = 3;
const int nx = 5;
const int ny = 5;
const int nz = 5;
// Create the Silo file
DBfile* file = DBCreate("sample.silo", DB_CLOBBER, DB_LOCAL, NULL, DB_PDB);
// names of the coordinates
char* coordnames[dim] = {"x", "y", "z"};
// dimensions
int dimensions[dim];
dimensions[0] = nx;
dimensions[1] = ny;
dimensions[2] = nz;
// the coordinate arrays
float node_x[nx] = {-2.0, -1.0, 0.0, 1.0, 2.0};
float node_y[ny] = {-2.0, -1.0, 0.0, 1.0, 2.0};
float node_z[nz] = {-2.0, -1.0, 0.0, 1.0, 2.0};
// the array of coordinate arrays
float* coordinates[dim];
coordinates[0] = node_x;
coordinates[1] = node_y;
coordinates[2] = node_z;
float data[nx * ny * nz];
for(int k = 0; k < nz; ++k){
for(int j = 0; j < ny; ++j){
for(int i = 0; i < nx; ++i){
data[i + nx*j + nx*ny*k] = i + nx * j + nx * ny * k;
}
}
}
// boost::multi_array<float, 3> data(boost::extents[nx][ny][nz]);
float* vars[] = {data};
DBPutQuadmesh(file, "mesh1", coordnames, coordinates, dimensions, dim, DB_FLOAT, DB_COLLINEAR, NULL);
DBPutQuadvar(file, "var1", "mesh1", 1, varnames, vars, dimensions, dim, NULL, 0, DB_FLOAT, DB_NODECENT, NULL);
// Close the Silo file
DBClose(file);
return 0;
}
とすればOKです。
Fortran形式のメモリ配置(列指向, column-major)のデータを渡す場合
BLASなどの計算用ライブラリを用いるためにデータの格納順を列指向にしていた場合、
DBOPT_MAJORORDERに1を渡せばよい。
...
// make options list
DBoptlist* optList = DBMakeOptlist(2);
DBAddOption(optList, DBOPT_UNITS, (void*)"V"); // 単位を与える場合
DBAddOption(optList, DBOPT_MAJORORDER, (void*)1); // column-major (Fortran) order
// tdArrayは boost::multi_array<double, 3> tdArray(boost::extents[nx][ny][nz], boost::fortran_storage_order());
double* vars[] = {tdArray.data()};
DBPutQuadmesh(file, "mesh1", coordnames, coordinates, dimensions, dim, DB_FLOAT, DB_COLLINEAR, NULL);
DBPutQuadvar(file, "var1", "mesh1", 1, varnames, vars, dimensions, dim, NULL, 0, DB_DOUBLE, DB_NODECENT, optList);
...
以上。
AMRで出力されたネストされたグリッド上のデータなどの出力法についてもそのうち書きたい。