소개

격자 계산 프로그램 생성 언어 Formura 을 사용해 본다.

그 3에서는 2차원 열전도 방정식(확산 방정식)을 풀어 보았다. 확산 방정식까지 오면, 조금 수정하는 것만으로 반응 확산 방정식을 풀 수 있다. 즉시 시도해 보자.

출처는 다음과 같습니다.

기주 b. 코 m/카이초 256/fm st

그 1 설치 및 컴파일까지

그 2 1차원 열전도 방정식

그 3 2차원 열전도 방정식

그 4 반응 확산 방정식 (Gray-Scott 계열) ← 이마 코코

Gray-Scott 방정식

반응 확산 방정식에는 다양한 것이 있지만, 비교적 식이 간단하고 결과가 재미있는 Gray-Scott계를 사용한다. 그 방정식은 다음과 같다.

\begin{align} 
\frac{\partial u}{\partial t} &= D_u \Delta u - uv^2 + F(1-u) \\
\frac{\partial v}{\partial t} &= D_v \Delta v + uv^2 - (F+k)v
\end{align}

여기서 $F$나 $k$는 정수, $D_u,D_v$는 확산 계수이다. 이 시스템에는 $u$와 $v$의 두 종류가 있기 때문에 그것을 Formura에서 다룬다.

fmr 파일

YAML 파일을 수정할 필요가 없습니다.

우선, 늘어난 분의 정수를 추가해 두자.

gs.fmr

double :: dt = 0.2
double :: Du = 0.05
double :: Dv = 0.1 # add
double :: F = 0.04 # add
double :: k = 0.06076 #add

쓰는 것을 잊었지만, # 이후는 코멘트로서 다루어진다.

초기화 함수 등으로, 2개의 상태 변수를 취급하는 경우는 튜플을 사용한다. 예를 들어, 초기화 함수는 다음과 같습니다.

gs.fmr

begin function (u,v) = init()
  double [] :: u,v
  u[i,j] = if isCenter(i,j,3) then 0.7 else 0.0
  v[i,j] = if isCenter(i,j,6) then 0.9 else 0.0
end function

확산 부분은 그대로. 역학 시스템의 부분은 다음과 같이 정의된다.

gs.fmr

calcU = fun(u,v) (u*u*v - (F+k)*u)
calcV = fun(u,v) (-u*u*v + F*(1.0-v))

각각 반응 확산 방정식의 '반응' 부분, 즉 오른쪽의 식을 그대로 Formura에 떨어뜨렸을 뿐이다.

시간 발전 부분을 쓰는 것도 어렵지 않다고 생각한다.

gs.fmr

begin function (u2,v2) = step(u,v)
  du = Du*diff(u)
  dv = Dv*diff(v)
  du = du + calcU(u,v)
  dv = dv + calcV(u,v)
  u2 = u + du * dt
  v2 = v + dv * dt
end function

두 개의 상태 변수가있는 경우 튜플을 사용하고 중간 변수를 사용할 수 있다는 것을 알면 위의 내용을 이해하는 것이 쉽습니다.

전체를 정리하면 이런 느낌이 된다.

gs.fmr

dimension :: 2
axes :: x,y

double :: dt = 0.2
double :: Du = 0.05
double :: Dv = 0.1 # add
double :: F = 0.04 # add
double :: k = 0.06076 #add

extern function :: fabs

isCenter = fun(i,j,w) (fabs(total_grid_x/2-i) < w) && (fabs(total_grid_y/2-j) < w)

diff = fun(q) (q[i+1,j] + q[i-1,j] + q[i,j+1] + q[i,j-1] - 4.0*q[i,j])

begin function (u,v) = init()
  double [] :: u,v
  u[i,j] = if isCenter(i,j,3) then 0.7 else 0.0
  v[i,j] = if isCenter(i,j,6) then 0.9 else 0.0
end function

calcU = fun(u,v) (u*u*v - (F+k)*u)
calcV = fun(u,v) (-u*u*v + F*(1.0-v))

begin function (u2,v2) = step(u,v)
  du = Du*diff(u)
  dv = Dv*diff(v)
  du = du + calcU(u,v)
  dv = dv + calcV(u,v)
  u2 = u + du * dt
  v2 = v + dv * dt
end function

「반응 확산 방정식을 차분화해 풀기」라고 하는 경우의 최저한의 기술이 되어 있는 것을 알 수 있을까 생각한다.

main 함수

main.cpp 의 수정은 불필요합니다만, 루프를 10000회로 해, 100회에 한번 덤프하도록 수정하자.

main.cpp

int main(int argc, char **argv) {
  Formura_Navi n;
  Formura_Init(&argc, &argv, &n);
  for (int i = 0; i < 10000; i++) {
    Formura_Forward(&n);
    if (i % 100 == 0) {
      dump(n);
    }
  }
  Formura_Finalize();
}

실행

그리고는 그대로 실행할 뿐이다.

$ formura gs.fmr
$ g++  main.cpp gs.c
$ rm -rf data
$ mkdir data
$ ./a.out
data/000.dat
data/001.dat
(snip)
data/098.dat
data/099.dat

플롯용 gnuplot 파일도 게재해 두자.

set term png
set xra [0:63]
set yra [0:63]
set view map
set size square
unset key

set cbrange[0:0.4]
do for[i=0:99:1]{
  input = sprintf("data/%03d.dat",i)
  output = sprintf("data/%03d.png",i)
  print output
  set out output
  sp input w pm3d
}

밝기 조정을 위해 cbrange를 수정했을 뿐이다. 가시화 결과는 이런 느낌이 된다.






애니메이션 GIF로 하면 이런 느낌.

convert -delay 10 -loop 0 *.png test.gif

결론

최저한의 기술로, 규칙 격자에 있어서의 차분법 코드를 출력해 주는 언어이며 프레임워크이기도 한 「Formura」를 사용해 보았다. 단순히 시뮬레이션 코드를 토하는 것 뿐만 아니라, OpenMP나 MPI에 의한 병렬화에도 대응하고 있어, 「교」풀 노드에서의 계산에도 성공, Gordon Bell상의 파이널리스트에도 선택되는 등하고 있지만, 병렬 화에 대해서는 다루지 않기 때문에 각자 시험되고 싶다.

Formura는 원래 무라주 쇼유키씨가 중심이 되어 개발한 것이다. 무라주씨는 「대단한 Haskell 즐겁게 배우자」의 번역자이기도 하며, PFN의 창립 멤버의 한 사람이기도 하고, ICPC(국제적인 경기 프로그래밍)에서도 활약하는 등, 업계에서는 널리 알려진 향후가 기대되는 젊은 재능이었다. 그 사카이가 정말 회개된다. 덧붙여 Formura는 계속해서 개발이 계속되고 있는 것 같다.