7年前に掲載したエントリー
myzyy.hatenablog.comに興味を持っていただいた方のご要望により、このときに使ったシミュレーションのプログラムを公開する。Fortran90で書き、gnuplotで可視化した。LinuxやMacOSなど、UNIX系のOSであれば、Fortran90はフリーソフトウェアのgfortranを使うことができると思う。gnuplotもフリーソフトウェアでありすぐにインストールできる。Windows10でもgfortranやgnuplotは使えるようだ。一連のプログラムをコンパイルして実行ファイルにした後、実行する。
(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
