はじめに
格子計算プログラム生成言語Formuraを使ってみる。
その1では、とりあえずFormuraのインストールと、Formuraが吐いたソースを読み込むmain関数を書いた。本稿では、簡単なシミュレーション例として一次元熱伝導方程式を解いてみる。前回と同じyamlとfmrを使い回すので、まだ読んでいない人は先にその1を参照してほしい。
初期条件の作成
初期条件として、中央に非ゼロの部分、あとは全部ゼロ、といったステップ関数的な温度を与える。そのために、gs.fmrに、「中央かどうか？」を判定する関数isCenterを定義しよう。
extern function :: fabs
isCenter = fun(i,w) (fabs(total_grid_x/2-i) < w)
Formuraは、外部関数としてC言語の関数を使える。ここではfabsを使っている。関数は、ラムダ式のようなイメージで定義できる。これは位置のインデックスiと、幅wを受け取り、中央から±wの幅だけTrueを返すような関数だ。
ここで、total_grid_xは自動生成された予約語で、x軸の全体のグリッド数である。
この関数を使って、「中央だけ非ゼロの値を入れる」という初期化関数は以下のように書ける。
begin function u = init()
double [] :: u
u[i] = if isCenter(i,3) then 0.7 else 0.0
end function
u[i]について、isCenter(i,3)がTrueなら0.7を、そうでなければ0.0を代入する式だ。これは添字iすべてについて実行される。
熱伝導方程式
Formuraファイル
次に、熱伝導方程式(拡散方程式)を定義しよう。この時、拡散係数Duや時間刻みdtを定義したくなるであろう。Formuraはこのようなスカラー変数も宣言できる。
double :: dt = 0.2
double :: Du = 0.05
熱伝導方程式の中心部分、ラプラシアンの計算を定義しよう。中央差分を取るだけなので簡単である。
diff = fun(q) (q[i+1] + q[i-1] - 2.0*q[i])
ラプラシアンが定義できれば、時間発展は以下のように簡単にかける。
begin function u2 = step(u)
u2 = u + diff(u) * dt
end function
これが熱伝導方程式を表していることは容易に見て取れるであろう。関数stepは、現在の状態変数ベクトルを受け取り、次の状態変数ベクトルを返す関数として定義する。Fortranの行列記述のように、添字を省略して記述できるし、状態変数ベクトルをそのまま先程定義した関数に入れて呼ぶことができる。
以上をすべてまとめると、Formuraファイルgs.fmrは以下のようになる。
dimension :: 1
axes :: x
double :: dt = 0.2
double :: Du = 0.05
extern function :: fabs
isCenter = fun(i,w) (fabs(total_grid_x/2-i) < w)
diff = fun(q) (q[i+1] + q[i-1] - 2.0*q[i])
begin function u = init()
double [] :: u
u[i] = if isCenter(i,3) then 0.7 else 0.0
end function
begin function u2 = step(u)
u2 = u + diff(u) * dt
end function
main関数
次に、main関数で時間発展を記述する。といってもFormura_Forward(&n)という関数を呼ぶだけで一回タイプステップが進むので簡単である。
毎ステップdumpという状態をダンプする関数を呼ぶなら、
int main(int argc, char **argv) {
Formura_Navi n;
Formura_Init(&argc, &argv, &n);
for (int i = 0; i < 100; i++) {
Formura_Forward(&n);
dump(n);
}
Formura_Finalize();
}
とかける。
状態の保存
Formuraでは、状態の保存には注意が必要だ。
まず、ベクトルのサイズはn.total_grid_xで取得できる。また、状態ベクトルuはformura_data.u[i]として取得できる。以上から、呼ばれるたびに状態ベクトルを保存する関数dumpを書くと、
void dump(Formura_Navi &n) {
char filename[256];
static int index = 0;
sprintf(filename, "data/%03d.dat", index);
index++;
std::cout << filename << std::endl;
std::ofstream ofs(filename);
for (int i = 0; i < n.total_grid_x; i++) {
double v = formura_data.u[i];
ofs << i << " ";
ofs << v << std::endl;
}
}
と書ける。しかし、このまま実行すると始点がずれてしまう。以下が実行結果だ。
$ formura gs.fmr
$ g++ main.cpp gs.c
$ mkdir data
$ ./a.out
data/000.dat
data/001.dat
(snip)
data/098.dat
data/099.dat
計算は正しくできていそうだが、場所がずれていることがわかるだろう。
これは実は、Formuraが毎回オフセットを変えて計算していることによる。このオフセットはn.offset_xで取得できるので、インデックスを以下のように補正してやれば良い。
int i2 = (i - n.offset_x + n.total_grid_x) % n.total_grid_x;
double v = formura_data.u[i2];
以下が実行結果だ。
軸が正しく修正できたことがわかるであろう。
以上をまとめたmain.cppは以下のようになる。
#include "gs.h"
#include <fstream>
#include <iostream>
void dump(Formura_Navi &n) {
char filename[256];
static int index = 0;
sprintf(filename, "data/%03d.dat", index);
index++;
std::cout << filename << std::endl;
std::ofstream ofs(filename);
for (int i = 0; i < n.total_grid_x; i++) {
int i2 = (i - n.offset_x + n.total_grid_x) % n.total_grid_x;
double v = formura_data.u[i2];
ofs << i << " ";
ofs << v << std::endl;
}
}
int main(int argc, char **argv) {
Formura_Navi n;
Formura_Init(&argc, &argv, &n);
for (int i = 0; i < 100; i++) {
Formura_Forward(&n);
dump(n);
}
Formura_Finalize();
}
まとめ
Formuraを使って、一次元熱伝導方程式を書いてみた。一度書いてしまえばmain.cppは使い回しができるので、後はFormuraファイル*.fmrを適当に修正するだけで自由に方程式を入れ替えることができる。
その3へ続く。