科学計算データの出力ライブラリSiloの利用法

# Siloとは
LLNLが開発している科学計算用のI/Oライブラリで、
様々なデータ形式用の出力に対応しています。
https://wci.llnl.gov/simulation/computer-codes/silo
同様にLLNLが開発している可視化ソフトウェアViSItがSilo形式を完全サポートしているのが特徴で、
データ出力から可視化までスムーズに行えます。
日本語情報が見つからなかったので、ビルドからC++での利用まで一通りまとめた内容を書いておきます。

# ビルド
```bash
$ 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章に示してあるサンプル通りに書いて出力してみる。

```C
#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する。
<img width="1288" alt="スクリーンショット 2016-12-21 16.47.47.png" src="https://qiita-image-store.s3.amazonaws.com/0/40906/a6ef1178-b6c2-be1a-a11d-703a2698da54.png">

meshが描画された。便利！！！

# 三次元の場合

DBPutQuadmeshの定義は

```C
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。

```C
#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で出力した図
<img width="916" alt="スクリーンショット 2016-12-21 17.35.33.png" src="https://qiita-image-store.s3.amazonaws.com/0/40906/2908b1b5-7997-73da-9f46-5c5e01dda7cc.png">

# 三次元データの出力
メッシュは出力できたので、対応するデータを出力します。
`Quadmesh`に対応するデータを入力するのは`DBPutQuadvar()`で、
その定義は

```C
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 length `nvars` containing pointers to arrays defining the values associated with each subvariable. For true edge- or face-centering (as opposed to `DB_EDGECENT` centering when `ndims` is 1 and `DB_FACECENT` centering when `ndims` is 2), each pointer here should point to an array that holds `ndims` sub-arrays, one for each of the i-, j-, k-oriented edges or i-, j-, k-intercepting faces, respectively. Read the description for more details.

サンプル

```C
#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を渡せばよい。

```C
...
        // 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で出力されたネストされたグリッド上のデータなどの出力法についてもそのうち書きたい。