7年前に掲載したエントリー
myzyy.hatenablog.comに興味を持っていただいた方のご要望により、このときに使ったシミュレーションのプログラムを公開する。Fortran90で書き、gnuplotで可視化した。LinuxやMacOSなど、UNIX系のOSであれば、Fortran90はフリーソフトウェアのgfortranを使うことができると思う。gnuplotもフリーソフトウェアでありすぐにインストールできる。Windows10でもgfortranやgnuplotは使えるようだ。一連のプログラムをコンパイルして実行ファイルにした後、実行する。
1) コンパイルのためのシェルスクリプト
(make_1dkl1.sh という名前でコピーペーストして保存していただきたい)
コマンドラインで
$ sh make_1dkl1.sh
と実行する。edge_city.1d.x という実行ファイルが作成される。
make_1dkl1.sh (1次元で求心力の及ぶ範囲が広い場合)
#! /bin/sh
PGM=edge_city.1d.x
F90=gfortran
F90OPT='-O3 -DKL1'
$F90 $F90OPT -c cal_inidens.F90
$F90 $F90OPT -c cal_potential.F90
$F90 $F90OPT -c edge_city.F90
$F90 $F90OPT -o $PGM *.o
rm -f *.o
make_1dkl2.sh (1次元で求心力の及ぶ範囲が狭い場合)
#! /bin/sh
PGM=edge_city.1d.x
F90=gfortran
F90OPT='-O3 -DKL2'
$F90 $F90OPT -c cal_inidens.F90
$F90 $F90OPT -c cal_potential.F90
$F90 $F90OPT -c edge_city.F90
$F90 $F90OPT -o $PGM *.o
rm -f *.o
make_2dkl1.sh (2次元で求心力の及ぶ範囲が広い場合)
#! /bin/sh
PGM=edge_city.2d.x
F90=gfortran
F90OPT='-O3 -DTWODIM -DKL1'
$F90 $F90OPT -c cal_inidens.F90
$F90 $F90OPT -c cal_potential.F90
$F90 $F90OPT -c edge_city.F90
$F90 $F90OPT -o $PGM *.o
rm -f *.o
他の場合も、プログラム名$PGMとプリプロセッサ―オプション-Dの場合分けでプログラムを別々に作成する。
2) プログラム
a) edge_city.F90 メインプログラム
ルンゲ・クッタ法で微分方程式の時間変化を離散化して表現(差分化)している。
プリプロセッサ―オプションでケース分けをしている。
-DKL1 (求心力の及ぶ反映が広い場合)
-DKL2 (求心力の及ぶ範囲が狭い場合)
-DTWODIM (2次元の場合)
-DKL2A47 (KL2の場合で、中心に求心力がはたらきその大きさがA=4.7)
-DKL2A48 (KL2の場合で、中心に求心力がはたらきその大きさがA=4.8)
program edge_city
!
#ifdef TWODIM
integer, parameter :: im = 24, jm = 24
#else
integer, parameter :: im = 24, jm = 1 ! 1D; also see cal_potential.f90
#endif
!
real, parameter :: dt = 0.1
real, parameter :: pi = 3.14159265
!
real :: smooth, scale, presum, nowsum
real :: A(im,jm), B(im,jm), dens(im,jm), r1, r2, gamma
real :: potential(im,jm), potential_mean
real :: tend1(im,jm), tend2(im,jm), tend3(im,jm), tend4(im,jm)
real :: dens1(im,jm), dens2(im,jm), dens3(im,jm), dens4(im,jm)
real :: mask(im,jm), dist, radius
integer :: step, stepend
character (len=5) :: cha
!
mask(:,:) = 1.
radius = im/2
do j = 1, jm
do i = 1, im
dist = sqrt(float((i-im/2)**2+(j-jm/2)**2))
if(dist > radius) mask(i,j) = 0.
end do
end do
!
stepend = 100000
scale = 1.
smooth = 5.
gamma = 1.
A = 0.2
B = 1.0
r1 = 1.0
r2 = 1.0
#ifdef KL1
r1 = 1.4
r2 = 0.2
#endif
#ifdef KL2
r1 = 2.8
r2 = 0.4
#endif
#ifdef KL2A48
write(6,*) 'center A KL2C4'
do j = 1, jm
do i = 1, im
A(i,j) = 4.8 * exp(-sqrt(float((i-im/2)**2+(j-jm/2)**2))*2.*pi/im*2.8)
end do
end do
#endif
#ifdef KL2A47
write(6,*) 'center A KL2C5'
do j = 1, jm
do i = 1, im
A(i,j) = 4.7 * exp(-sqrt(float((i-im/2)**2+(j-jm/2)**2))*2.*pi/im*2.8)
end do
end do
#endif
do j = jm, 1, -1
write(6,'(24(x,i3))') ((int(A(i,j)*10.)),i=1,im,1)
end do
write(6,*) 'r1 r2 A B ', r1, r2, A(1,1), B(1,1)
write(6,*) 'r1/r2 A/B r2/r1 ', r1/r2, A(1,1)/B(1,1), r2/r1
call cal_inidens(dens,mask,smooth,scale,im,jm)
do j = jm, 1, -1
write(6,'(24(x,i3))') ((int(dens(i,j)*100.)),i=1,im,1)
end do
!
presum = 10.
do step = 1, stepend
call cal_potential(potential,potential_mean,dens,mask,A,B,r1,r2,im,jm)
tend1(:,:) = gamma * (potential(:,:) - potential_mean) * dens(:,:)
dens2(:,:) = dens(:,:) + 0.5 * dt * tend1(:,:)
call cal_potential(potential,potential_mean,dens2,mask,A,B,r1,r2,im,jm)
tend2(:,:) = gamma * (potential(:,:) - potential_mean) * dens2(:,:)
dens3(:,:) = dens(:,:) + 0.5 * dt * tend2(:,:)
call cal_potential(potential,potential_mean,dens3,mask,A,B,r1,r2,im,jm)
tend3(:,:) = gamma * (potential(:,:) - potential_mean) * dens3(:,:)
dens4(:,:) = dens(:,:) + dt * tend3(:,:)
call cal_potential(potential,potential_mean,dens4,mask,A,B,r1,r2,im,jm)
tend4(:,:) = gamma * (potential(:,:) - potential_mean) * dens4(:,:)
dens(:,:) = dens(:,:) + dt/6. * &
& (tend1(:,:) + 2.*tend2(:,:) + 2.*tend3(:,:) + tend4(:,:))
dens(:,:) = dens(:,:)/sum(dens)
nowsum = 0.5*(sum(dt * abs(tend1)) + presum)
if(nowsum < 1.e-6) exit
if(mod(step,100) == 0) then
write(6,*) 'step ', step, nowsum, presum, sum(dens)
do j = jm, 1, -1
write(6,'(24(x,i3))') ((int(dens(i,j)*100.)),i=1,im,1)
end do
call flush(6)
end if
presum = nowsum
if(mod(step,100) == 0) then
write(cha,'(i5.5)') int(step*dt)
open(10,file='stepdt.'//cha//'.txt',form='formatted')
do j = jm, 1, -1
write(10,'(24(x,f10.7))') (dens(i,j),i=1,im,1)
end do
close(10)
end if
end do
!
write(6,*) '### dt end ', sum(dt * abs(tend1))
write(6,*) '### step tend dens ', step, sum(abs(tend4)), sum(dens)
do j = jm, 1, -1
write(6,'(24(x,i3))') ((int(dens(i,j)*100.)),i=1,im,1)
end do
open(10,file='stepdt.final.txt',form='formatted')
do j = jm, 1, -1
write(10,'(24(x,f10.7))') (dens(i,j),i=1,im,1)
end do
close(10)
!
end
b) cal_inidens.F90 初期値作成
subroutine cal_inidens(dens,mask,smooth,scale,im,jm)
!
integer, intent(in) :: im, jm
real, intent(in) :: smooth, scale, mask(im,jm)
real, intent(out) :: dens(im,jm)
!
real :: ran(im,jm), sum
integer :: i, j
!
sum = 0.
call random_number(ran)
do j = 1, jm
do i = 1, im
dens(i,j) = (ran(i,j) + smooth) * scale * mask(i,j)
sum = sum + dens(i,j)
end do
end do
dens(:,:) = dens(:,:) / sum
!
return
!
end subroutine cal_inidens
c) cal_potential.F90 動学モデル(1)の右辺にある項(2)の計算を行うサブルーチン
subroutine cal_potential(potential,potential_mean,dens,mask,A,B,r1,r2,im,jm)
!
integer, intent(in) :: im, jm
real, intent(in) :: A(im,jm), B(im,jm), dens(im,jm), mask(im,jm), r1, r2
real, intent(out) :: potential(im,jm), potential_mean
!
real :: dist, dist1, dist2, dist3
integer i, j, m, n
real, parameter :: pi = 3.14159265
!
potential(:,:) = 0.
potential_mean = 0.
do j = 1, jm
do i = 1, im
do n = 1, jm
do m = 1, im
#ifdef TWODIM
dist = sqrt(float((n - j)**2 + (m - i)**2))*2.*pi/im
#else
dist1 = sqrt(float((n - j)**2 + (m - i)**2))*2.*pi/im
dist2 = sqrt(float((n - j)**2 + (m-im -i)**2))*2.*pi/im
dist3 = sqrt(float((n - j)**2 + (m+im -i)**2))*2.*pi/im
dist = min(dist1,dist2,dist3)
#endif
potential(i,j) = potential(i,j) &
& +(A(i,j) * exp(-r1 * dist) &
& - B(i,j) * exp(-r2 * dist)) &
& * dens(m,n) * mask(i,j)
end do
end do
potential(i,j) = mask(i,j) * potential(i,j)
potential_mean = potential_mean + potential(i,j) * dens(i,j)
end do
end do
!
return
!
end subroutine cal_potential
3) 可視化シェルスクリプト
a) tplot.sh (1次元) ($ sh tplot.sh を実行すると、時間変化を含めてすべての状態を可視化する)
#! /bin/sh
cat stepdt.0*.txt > stepdt.all.txt
EPS=tplot.eps
PNG=tplot.png
STEP=all
GNP=tplot.gnp
cat <$GNP
set output '$EPS'
set term postscript eps enhanced color 12
set view 60,300
set ticslevel 2
set pm3d at sb
set cntrparam bspline
set palette rgbformulae 22, 13, -31
set xtics font 'Helvetica,24'
set xlabel 'position'
set xlabel font 'Helvetica,24'
set ytics font 'Helvetica,24'
set ylabel 'iteration (X 100)'
set ylabel font 'Helvetica,24'
set nokey
splot 'stepdt.$STEP.txt' matrix with lines
EOF
gnuplot $GNP
convert -trim $EPS $PNG
b) splot.sh (2次元) ($ sh splot.sh 002760 などとして実行する。各ステップ数のスナップショットを描画する)
#! /bin/sh
STEP=$1
EPS=splot_$STEP.eps
PNG=splot_$STEP.png
GNP=splot.gnp
cat <$GNP
set output '$EPS'
set term postscript eps color
set view 60,15
set ticslevel 2
set pm3d at sb
set cntrparam bspline
set palette rgbformulae 22, 13, -31
set xlabel 'position'
set xlabel font 'Helvetica,24'
set xtics font 'Helvetica,24'
set ylabel 'position'
set ylabel font 'Helvetica,24'
set ytics font 'Helvetica,24'
set nokey
splot 'stepdt.$STEP.txt' matrix with lines
EOF
gnuplot $GNP
convert -trim $EPS $PNG