Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
Help us understand the problem. What is going on with this article?

格子計算プログラム生成言語Formuraを使ってみる その2

More than 1 year has passed since last update.

はじめに

格子計算プログラム生成言語Formuraを使ってみる。

その1では、とりあえずFormuraのインストールと、Formuraが吐いたソースを読み込むmain関数を書いた。本稿では、簡単なシミュレーション例として一次元熱伝導方程式を解いてみる。前回と同じyamlとfmrを使い回すので、まだ読んでいない人は先にその1を参照してほしい。

  • その1 インストールとコンパイルまで
  • その2 一次元熱伝導方程式 ←イマココ
  • その3 二次元熱伝導方程式
  • その4 反応拡散方程式(Gray-Scott系)

初期条件の作成

初期条件として、中央に非ゼロの部分、あとは全部ゼロ、といったステップ関数的な温度を与える。そのために、gs.fmrに、「中央かどうか?」を判定する関数isCenterを定義しよう。

gs.fmr
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軸の全体のグリッド数である。

この関数を使って、「中央だけ非ゼロの値を入れる」という初期化関数は以下のように書ける。

gs.fmr
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はこのようなスカラー変数も宣言できる。

gs.fmr
double :: dt = 0.2
double :: Du = 0.05

熱伝導方程式の中心部分、ラプラシアンの計算を定義しよう。中央差分を取るだけなので簡単である。

gs.fmr
diff = fun(q) (q[i+1] + q[i-1] - 2.0*q[i])

ラプラシアンが定義できれば、時間発展は以下のように簡単にかける。

gs.fmr
begin function u2 = step(u)
  u2 = u + diff(u) * dt
end function

これが熱伝導方程式を表していることは容易に見て取れるであろう。関数stepは、現在の状態変数ベクトルを受け取り、次の状態変数ベクトルを返す関数として定義する。Fortranの行列記述のように、添字を省略して記述できるし、状態変数ベクトルをそのまま先程定義した関数に入れて呼ぶことができる。

以上をすべてまとめると、Formuraファイルgs.fmrは以下のようになる。

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で取得できる。また、状態ベクトルuformura_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

image.png

計算は正しくできていそうだが、場所がずれていることがわかるだろう。

これは実は、Formuraが毎回オフセットを変えて計算していることによる。このオフセットはn.offset_xで取得できるので、インデックスを以下のように補正してやれば良い。

    int i2 = (i - n.offset_x + n.total_grid_x) % n.total_grid_x;
    double v = formura_data.u[i2];

以下が実行結果だ。

image.png

軸が正しく修正できたことがわかるであろう。

以上をまとめたmain.cppは以下のようになる。

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へ続く。

kaityo256
記事中に明示されていない場合、私の記事はCC-BY 4.0で、記事に含まれるソースコードはMITライセンスで公開します。
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away