FortranのコードをJuliaへ移植してみる
マンデルブロ集合
こちらの記事にはさまざまなFortranのバージョンで書かれたマンデルブロ集合のコードが紹介されています。
Fortran90
この記事の中のFortran90のコードを引用します。コードは
PROGRAM mandel
IMPLICIT NONE
INTEGER, PARAMETER :: nx = 61, ny = 31, maxiter = 90
REAL , PARAMETER :: x0 = -2.0, x1 = 2.0, y0 = -2.0, y1 = 2.0
INTEGER :: ix, iy, iter, mandelbrot(nx, ny)
REAL :: x, y
COMPLEX :: z, c
DO iy = 1, ny
y = y0 + (y1 - y0) * (iy - 1) / REAL(ny - 1)
DO ix = 1, nx
x = x0 + (x1 - x0) * (ix - 1) / REAL(nx - 1)
c = CMPLX(x, y)
z = (0.0, 0.0)
DO iter = 0, maxiter
z = z * z + c
IF (ABS(z) > 2.0) EXIT
END DO
mandelbrot(ix, iy) = iter
END DO
END DO
!
DO iy = 1, ny
PRINT '(61I1)', (mandelbrot(:, iy) + 9) / 10
END DO
STOP
END PROGRAM mandelです。このコードをgfortranでコンパイルし、実行しますと、
0000000000000000000000000000001000000000000000000000000000000
0000000000000000000011111111111111111111100000000000000000000
0000000000000000111111111111111111111111111110000000000000000
0000000000001111111111111111111111111111111111111000000000000
0000000000111111111111111111111111111111111111111110000000000
0000000011111111111111111111111111111111111111111111100000000
0000001111111111111111111111111111111111111111111111111000000
0000011111111111111111111111111111111111111111111111111100000
0000111111111111111111111111211111111111111111111111111110000
000111111111111111111111111***2111111111111111111111111111000
0011111111111111111113312223*72671111111111111111111111111100
0011111111111111111111*************61111111111111111111111100
011111111111111111112***************1111111111111111111111110
0111111111112*4*3412****************2111111111111111111111110
011111111212************************5111111111111111111111110
**********************************311111111111111111111111111
011111111212************************5111111111111111111111110
0111111111112*4*3412****************2111111111111111111111110
011111111111111111112***************1111111111111111111111110
0011111111111111111111*************61111111111111111111111100
0011111111111111111113312223*72661111111111111111111111111100
000111111111111111111111111***2111111111111111111111111111000
0000111111111111111111111111211111111111111111111111111110000
0000011111111111111111111111111111111111111111111111111100000
0000001111111111111111111111111111111111111111111111111000000
0000000011111111111111111111111111111111111111111111100000000
0000000000111111111111111111111111111111111111111110000000000
0000000000001111111111111111111111111111111111111000000000000
0000000000000000111111111111111111111111111110000000000000000
0000000000000000000011111111111111111111100000000000000000000
0000000000000000000000000000001000000000000000000000000000000という出力が得られます。この出力を再現するようなJuliaのコードを書いてみましょう。
作ったJuliaのコードは
function mandel()
nx = 61; ny = 31; maxiter = 90
x0 = -2.0; x1 = 2.0; y0 = -2.0; y1 = 2.0
mandelbrot = zeros(Int64,nx, ny)
mandelbrot .= maxiter
for iy=1:ny
y = y0 + (y1 - y0) * (iy - 1) / real(ny - 1)
for ix=1:nx
x = x0 + (x1 - x0) * (ix - 1) / real(nx - 1)
c = x + im*y
z = 0im
for iter=0:maxiter
z = z * z + c
if abs(z) > 2
mandelbrot[ix,iy] = iter
break
end
end
end
end
for iy=1:ny
for ix=1:nx
if mandelbrot[ix,iy] == maxiter
print("*")
else
print((mandelbrot[ix,iy]+9) ÷ 10)
end
end
println("\t")
end
end
mandel()です。
ポイントは
program mandelをfunction mandel()に変更;を使い一行に複数の文を入れた- 配列
mandelbrotをmandelbrot = zeros(Int64,nx, ny)で定義 - Fortranでループを抜けるexitは、対応するbreakに変更
- Fortranは
do iter = 0, maxiterのループの途中で抜けた場合iterには抜けた時の値が入っているが、Juliaではiterのスコープがローカルなので外では未定義になっている。そのため、mandelbrot[ix,iy] = iterとした iterのループが最大まで到達した時に、Fortranではmandelbrot(ix, iy)の中身はループの最後の値maxiterになるが、Juliaでは前述のようにiterはローカルスコープの変数なので値が入らない。なので、最初にmandelbrot .= maxiterと初期化した
などです。見比べるとはっきりとわかりますが、Fortran90のコードとJuliaのコードはよく似ています。
Fortran90 モジュールを使用
次は、モジュールを使用したコードをJuliaコードに変化させてみます。 コードはこちらの記事からの引用しまして、
MODULE m_mandel
IMPLICIT NONE
INTEGER, PARAMETER :: maxiter = 90
CONTAINS
INTEGER FUNCTION mandel(c)
COMPLEX, INTENT(IN) :: c
COMPLEX :: z
z = (0.0, 0.0)
DO mandel = 0, maxiter
z = z * z + c
IF (ABS(z) > 2.0) EXIT
END DO
END FUNCTION mandel
END MODULE m_mandel
PROGRAM mandel_main
USE m_mandel
IMPLICIT NONE
INTEGER, PARAMETER :: nx = 61, ny = 31
REAL , PARAMETER :: x0 = -2.0, x1 = 2.0, y0 = -2.0, y1 = 2.0
INTEGER :: ix, iy, iter, mandelbrot(nx, ny)
REAL :: x, y
DO iy = 1, ny
y = y0 + (y1 - y0) * (iy - 1) / REAL(ny - 1)
DO ix = 1, nx
x = x0 + (x1 - x0) * (ix - 1) / REAL(nx - 1)
mandelbrot(ix, iy) = mandel(CMPLX(x, y))
END DO
END DO
!
DO iy = 1, ny
PRINT '(61i1)', (mandelbrot(:, iy) + 9) / 10
END DO
STOP
END PROGRAM mandel_mainというものです。
このコードをJuliaに移植してみますと、
module M_mandel
const maxiter = 90
export mandel,maxiter
function mandel(c)
z = 0im
count = maxiter
for i=0:maxiter
z = z * z + c
if abs(z) > 2
count = i
break
end
end
return count
end
end
using .M_mandel
function mandel_main()
nx = 61;ny = 31
x0 = -2.0;x1 = 2.0;y0 = -2.0;y1 = 2.0
mandelbrot = zeros(Int64,nx, ny)
for iy=1:ny
y = y0 + (y1 - y0) * (iy - 1) / real(ny - 1)
for ix=1:nx
x = x0 + (x1 - x0) * (ix - 1) / real(nx - 1)
mandelbrot[ix, iy] = mandel(x+im*y)
end
end
for iy=1:ny
for ix=1:nx
if mandelbrot[ix,iy] == maxiter
print("*")
else
print((mandelbrot[ix,iy]+9) ÷ 10)
end
end
println("\t")
end
end
mandel_main()となります。
一つ目のFortran90コードで気をつけたポイントの他には、
MODULEはmoduleに。モジュールの名前はJuliaでは最初の一文字を大文字にする慣例があるのでそれにならったUSE m_mandelはusing .M_mandelに。関数の中ではなく、外に定義する- Fortranではモジュール内で宣言した変数や関数はデフォルトで
share属性を持っておりuseですぐに使えるが、Juliaの場合に同様にしたい場合にはexportをモジュール内で使っておく。使わない場合はM_mandel.mandel(x+im*y)とアクセスすることになる
ということに気をつけました。
Fortran2003
次は、同じ記事のFortran2003のコードをJuliaに移植してみます。Fortranのコードを引用すると
program mandel
implicit none
integer, parameter :: maxiter = 90, nx = 61, ny = 31
real , parameter :: x0 = -2.0, x1 = 2.0, y0 = -2.0, y1 = 2.0
character(len = 10), parameter :: text = '|+o0O.@*#-'
integer :: ix, iy, iter, mandelbrot(nx, ny) = 0, m(nx)
real, allocatable :: x(:), y(:)
complex :: c(nx, ny), z(nx, ny) = (0.0, 0.0)
x = [( (x1 - x0) / (nx - 1) * (ix - 1) + x0, ix = 1, nx )]
y = [( (y1 - y0) / (ny - 1) * (iy - 1) + y0, iy = 1, ny )]
forall (ix = 1:nx, iy = 1:ny) c(ix, iy) = cmplx(x(ix), y(iy))
do iter = 0, maxiter
where (abs(z) <= 2.0)
z = z * z + c
mandelbrot = mandelbrot + 1
end where
end do
do iy = 1, ny
m = (mandelbrot(:, iy) + 8) / 10 + 1
write(*, '(61a1)') (text(m(ix):m(ix)), ix = 1, nx)
end do
stop
end program mandelとなります。このコードをコンパイルして実行すると、
||||||||||||||||||||||||||||||+||||||||||||||||||||||||||||||
||||||||||||||||||||+++++++++++++++++++++||||||||||||||||||||
||||||||||||||||+++++++++++++++++++++++++++++||||||||||||||||
||||||||||||+++++++++++++++++++++++++++++++++++++||||||||||||
||||||||||+++++++++++++++++++++++++++++++++++++++++||||||||||
||||||||+++++++++++++++++++++++++++++++++++++++++++++||||||||
||||||+++++++++++++++++++++++++++++++++++++++++++++++++||||||
|||||+++++++++++++++++++++++++++++++++++++++++++++++++++|||||
||||++++++++++++++++++++++++o++++++++++++++++++++++++++++||||
|||++++++++++++++++++++++++---o+++++++++++++++++++++++++++|||
||+++++++++++++++++++00+ooo0-*o@*++++++++++++++++++++++++++||
||++++++++++++++++++++-------------@+++++++++++++++++++++++||
|+++++++++++++++++++o---------------++++++++++++++++++++++++|
|+++++++++++o-O-0O+o----------------o+++++++++++++++++++++++|
|++++++++o+o------------------------.+++++++++++++++++++++++|
----------------------------------0++++++++++++++++++++++++++
|++++++++o+o------------------------.+++++++++++++++++++++++|
|+++++++++++o-O-0O+o----------------o+++++++++++++++++++++++|
|+++++++++++++++++++o---------------++++++++++++++++++++++++|
||++++++++++++++++++++-------------@+++++++++++++++++++++++||
||+++++++++++++++++++00+ooo0-*o@@++++++++++++++++++++++++++||
|||++++++++++++++++++++++++---o+++++++++++++++++++++++++++|||
||||++++++++++++++++++++++++o++++++++++++++++++++++++++++||||
|||||+++++++++++++++++++++++++++++++++++++++++++++++++++|||||
||||||++++++++++++++++++++++++++++++++++++++++++++++++|||||||
||||||||+++++++++++++++++++++++++++++++++++++++++++++||||||||
||||||||||+++++++++++++++++++++++++++++++++++++++++||||||||||
||||||||||||++++++++++++++++++++++++++++++++++++|||||||||||||
||||||||||||||||+++++++++++++++++++++++++++++||||||||||||||||
||||||||||||||||||||+++++++++++++++++++++||||||||||||||||||||
||||||||||||||||||||||||||||||+||||||||||||||||||||||||||||||となります。
対応するJuliaのコードは
function mandel()
maxiter = 90;nx = 61; ny = 31
x0 = -2.0; x1 = 2.0; y0 = -2.0; y1 = 2.0
text = "|+o0O.@*#-"
z = zeros(ComplexF64,nx,ny)
mandelbrot = zeros(Int64,nx,ny)
x = [ (x1 - x0) / (nx - 1) * (ix - 1) + x0 for ix = 1:nx]
y = [ (y1 - y0) / (ny - 1) * (iy - 1) + y0 for iy = 1:ny]
c = [ x[ix]+im*y[iy] for ix=1:nx,iy=1:ny]
for iter=0:maxiter
for i=1:length(z)
if abs(z[i]) <= 2
z[i] = z[i]*z[i]+c[i]
mandelbrot[i] += 1
end
end
end
for iy=1:ny
m = (mandelbrot[:, iy] .+ 8) .÷ 10 .+ 1
[print(text[m[ix]:m[ix]]) for ix=1:nx]
println("\t")
end
end
mandel()です。Fortranでのwhere構文（配列の各要素に条件式を当てはめてループする構文）のJulia対応が分かりませんでしたので、仕方がないのでforループを使いました。ただし、そのまま書いても芸がありませんので、2次元配列のzなどを1次元配列として取り出す方法を使ってみました。Juliaでは2次元配列をメモリ格納順に並べた1次元のインデックスを用いて値を操作することができます。
今回は
- for文のリスト内包表記
[ (x1 - x0) / (nx - 1) * (ix - 1) + x0 for ix = 1:nx]の使用 - Fortranの
forallの代わりに二重forループのリスト内包表記の使用 - 配列の各要素にまとめて演算するために
.を使ったブロードキャストを使用（m = (mandelbrot[:, iy] .+ 8) .÷ 10 .+ 1）
を行いました。
Fortran2003 オブジェクト指向
次は、同じ記事のオブジェクト指向で書かれたFortranコードをJuliaに移植してみます。 コードは
module m_bmp
implicit none
type :: t_bmp_file_header
sequence ! 14bytes
character(len = 2) :: bfType = 'BM' !integer(2) :: bfType = transfer('BM', 0_2, 1) ! BitMap
integer(4) :: bfSize ! file size in bytes
integer(2) :: bfReserved1 = 0 ! always 0
integer(2) :: bfReserved2 = 0 ! always 0
integer(4) :: bfOffBits
end type t_bmp_file_header
!
type :: t_bmp_info_header
sequence
integer(4) :: biSize = Z'28' ! size of bmp_info_header ; 40bytes
integer(4) :: biWidth
integer(4) :: biHeight
integer(2) :: biPlanes = 1 ! always 1
integer(2) :: biBitCount
integer(4) :: biCompression = 0 !0:nocompression,1:8bitRLE,2:4bitRLE,3:bitfield
integer(4) :: biSizeImage
integer(4) :: biXPelsPerMeter = 3780 ! 96dpi
integer(4) :: biYPelsPerMeter = 3780 ! 96dpi
integer(4) :: biClrUsed = 0
integer(4) :: biClrImportant = 0
end type t_bmp_info_header
!
type :: t_rgb
sequence
character :: b, g, r ! order is b g r
end type t_rgb
!
type :: t_bmp(nx, ny)
integer, len:: nx, ny
type(t_rgb) :: rgb(nx, ny)
contains
procedure :: wr => wr_bmp
procedure :: pr_bmp
generic :: write(formatted) => pr_bmp
end type
contains
subroutine wr_bmp(bmp, fn)
class(t_bmp(*, *)), intent(in) :: bmp
character(len = *), intent(in) :: fn
type(t_bmp_file_header) :: bmp_file_header
type(t_bmp_info_header) :: bmp_info_header
bmp_file_header%bfSize = 14 + bmp_info_header%biSize + 0 + bmp%nx * bmp%ny * 3
bmp_file_header%bfOffBits = 14 + bmp_info_header%biSize
bmp_info_header%biWidth = bmp%nx ! nx shouold be a multiple of 4
bmp_info_header%biHeight = bmp%ny
bmp_info_header%biBitCount = 24 ! color depth 24bits
bmp_info_header%biSizeImage = bmp%nx * bmp%ny * 3
open(9, file = fn//'.bmp', form = 'binary', status = 'unknown')
write(9) bmp_file_header
write(9) bmp_info_header
write(9) bmp%rgb
close(9)
return
end subroutine wr_bmp
! convert to t_RGB
pure elemental type(t_rgb) function to_rgb(ir, ig, ib)
integer, intent(in) :: ir, ig, ib
to_rgb = t_rgb(achar(ib), achar(ig), achar(ir))
end function to_rgb
subroutine pr_bmp(dtv, unit, iotype, vlist, io, iomsg)
class(t_bmp(*, *)), intent(in) :: dtv
integer, intent(in) :: unit
character(len = *), intent(in) :: iotype
integer, intent(in) :: vlist(:)
integer, intent(out) :: io
character(len = *), intent(in out) :: iomsg
character(len = 30) :: fmt
if (iotype == 'LISTDIRECTED') then
write(unit, *, iostat = io) 'nx =', dtv%nx, ', ny =', dtv%ny
end if
end subroutine pr_bmp
end module m_bmp
module m_mandel
implicit none
integer, parameter :: maxiter = 254
contains
pure elemental integer function mandel(c)
complex, intent(in) :: c
complex :: z
z = (0.0, 0.0)
do mandel = 0, maxiter
z = z * z + c
if (abs(z) > 2.0) exit
end do
end function mandel
end module m_mandel
program mandel_main
use m_bmp
use m_mandel
implicit none
integer, parameter :: nx = 640, ny = 640
real , parameter :: x0 = -2.0, x1 = 2.0, y0 = -2.0, y1 = 2.0
integer :: ix, iy, iter, mandelbrot(nx, ny)
real :: x(nx), y(ny)
complex :: c(nx, ny)
type(t_bmp(nx, ny)) :: bmp
!
forall (ix = 1:nx) x(ix) = x0 + (x1 - x0) * (ix - 1) / real(nx - 1)
forall (iy = 1:ny) y(iy) = y0 + (y1 - y0) * (iy - 1) / real(ny - 1)
forall (ix = 1:nx, iy = 1:ny) c(ix, iy) = cmplx(x(ix), y(iy))
mandelbrot = mandel(c)
!
bmp%rgb = to_rgb(255 - mandelbrot, 255 - mandelbrot, 255 - mandelbrot)
call bmp%wr('test')
print *, 'BMP size: ', bmp
stop
end program mandel_mainです。このコードですが、gfortranでコンパイルしようとしても失敗してしまいました。こちらによると、どうやら「Fortran2003 で導入されたパラメータ付派生型(parameterized derived type)」がgfortranがこの機能にフルに対応していないようです。ですので、このコードの実行を確認できていません。
追記
2020年12月からIntelのFortranコンパイラであるifortが無料で使えるようになっていたようです。こちらを参照してインストールしてみたところ、無事に上のコードをコンパイルできました。なお、Macでやる時にはこちらにあるように「 -L/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/lib」をつける必要があるかもしれません（LIBRARY_PATHにこれを追加しておく）。
記事によると、このコードは画像をBMP形式で出力してくれるようです。BMPのフォーマットは検索すると出てきましたので、その形式に従ってバイナリを出力する、というコードのようです。
これのJulia版は
module M_bmp
export to_rgb,T_bmp,wr_bmp
Base.@kwdef mutable struct T_bmp_file_header
bfType::String = "BM"
bfSize::Int32 = 0
bfReserved1::Int16 = 0
bfReserved2::Int16 = 0
bfOffBits::Int32 = 0
end
Base.@kwdef mutable struct T_bmp_info_header
biSize::Int32 = 0x28
biWidth::Int32 = 0
biHeight::Int32 = 0
biPlanes::Int16 = 1
biBitCount::Int16 = 0
biCompression::Int32 = 0 #0:nocompression,1:8bitRLE,2:4bitRLE,3:bitfield
biSizeImage::Int32 = 0
biXPelsPerMeter::Int32 = 3780 # 96dpi
biYPelsPerMeter::Int32 = 3780 # 96dpi
biClrUsed::Int32 = 0
biClrImportant::Int32 = 0
end
mutable struct T_rgb
b::UInt8
g::UInt8
r::UInt8
end
Base.@kwdef mutable struct T_bmp{nx,ny}
rgb::Array{T_rgb,2} = Array{T_rgb,2}(undef,nx,ny)
end
function Base.write(fp::IO,data::T_bmp_file_header)
write(fp,data.bfType)
write(fp,data.bfSize)
write(fp,data.bfReserved1)
write(fp,data.bfReserved2)
write(fp,data.bfOffBits)
end
function Base.write(fp::IO,data::T_bmp_info_header)
write(fp,data.biSize)
write(fp,data.biWidth)
write(fp,data.biHeight)
write(fp,data.biPlanes)
write(fp,data.biBitCount)
write(fp,data.biCompression)
write(fp,data.biSizeImage)
write(fp,data.biXPelsPerMeter)
write(fp,data.biYPelsPerMeter)
write(fp,data.biClrUsed)
write(fp,data.biClrImportant)
end
function Base.write(fp::IO,data::T_rgb)
write(fp,data.b)
write(fp,data.g)
write(fp,data.r)
end
function Base.write(fp::IO,data::T_bmp{nx,ny}) where {nx,ny}
for iy=1:ny
for ix=1:nx
#println(data.rgb[ix,iy])
write(fp,data.rgb[ix,iy])
end
end
end
function wr_bmp(bmp::T_bmp{nx,ny}, fn) where {nx,ny}
bmp_info_header = T_bmp_info_header()
bmp_file_header = T_bmp_file_header()
bmp_file_header.bfSize = 14 + bmp_info_header.biSize + 0 + nx * ny * 3
bmp_file_header.bfOffBits = 14 + bmp_info_header.biSize
bmp_info_header.biWidth = nx # nx shouold be a multiple of 4
bmp_info_header.biHeight = ny
bmp_info_header.biBitCount = 24 # color depth 24bits
bmp_info_header.biSizeImage = nx * ny * 3
fp = open(fn*".bmp","w")
write(fp,bmp_file_header)
write(fp,bmp_info_header)
write(fp,bmp)
end
function to_rgb(ir,ig,ib)
return T_rgb(htol(ir),htol(ig),htol(ib))
end
end
module M_mandel
const maxiter = 254
export mandel,maxiter
function mandel(c)
z = 0im
count = maxiter
for i=0:maxiter
z = z * z + c
if abs(z) > 2
count = i
break
end
end
return count
end
end
using .M_bmp
using .M_mandel
function mandel_main()
nx = 640;ny = 640
x0 = -2.0;x1 = 2.0; y0 = -2.0; y1 = 2.0
bmp = T_bmp{nx,ny}()
x = [ (x1 - x0) / (nx - 1) * (ix - 1) + x0 for ix = 1:nx]
y = [ (y1 - y0) / (ny - 1) * (iy - 1) + y0 for iy = 1:ny]
c = [ x[ix]+im*y[iy] for ix=1:nx,iy=1:ny]
mandelbrot = mandel.(c)
bmp.rgb = to_rgb.(255 .- mandelbrot,255 .- mandelbrot,255 .- mandelbrot)
wr_bmp(bmp,"test")
println("BMP size: $nx $ny")
end
mandel_main()となります。 実行した結果は得られる画像は
となります。
ポイントは
- Fortran2003のオブジェクト指向によるクラスはJuliaでは
mutable structに。また、クラスのメソッドはJuliaでは多重ディスパッチで定義。 - Fortranでは定義した構造体をそのままバイナリで書き出せたが、Juliaでの独自型は自分で書き出しを指定してやる必要があるため、標準のwriteを多重ディスパッチで機能拡張
- BMP形式の定義によると
t_rgbクラスに入れるべきは符号なし整数なので、JuliaではUInt8型（符号なし8ビット整数）に変更 mutable structでデフォルトの値を入れるためにBase.@kwdefを使用- Fortranでのパラメータ付派生型と同等の機能をJuliaではパラメトリック型で実装（
mutable struct T_bmp{nx,ny}）
です。オブジェクト指向なコードもすんなりJuliaに移植できました。
パッケージの利用
上で書いたJuliaのコードはFortranでBMP形式を扱うように直接バイナリを操作したものでした。Juliaの場合、自前で書かなくても関連するパッケージがあり、簡単に画像データを扱えます。ここではImagesというパッケージを使います。これをインストールするには、
julia -e 'using Pkg;Pkg.add("Images")'としてください。 これを使うと、
module M_mandel
const maxiter = 254
export mandel,maxiter
function mandel(c)
z = 0im
count = maxiter
for i=0:maxiter
z = z * z + c
if abs(z) > 2
count = i
break
end
end
return count
end
end
using .M_mandel
using Images
function mandel_main()
nx = 640;ny = 640
x0 = -2.0;x1 = 2.0; y0 = -2.0; y1 = 2.0
x = [ (x1 - x0) / (nx - 1) * (ix - 1) + x0 for ix = 1:nx]
y = [ (y1 - y0) / (ny - 1) * (iy - 1) + y0 for iy = 1:ny]
c = [ x[ix]+im*y[iy] for ix=1:nx,iy=1:ny]
mandelbrot = mandel.(c)
img = (255 .- mandelbrot)/255
save("test2.bmp",img')
save("test2.png",img')
println("BMP size: $nx $ny")
end
mandel_main()とすれば、BMP形式とPNG形式のグレースケールの画像を出力できます。
ここでの追加のポイントとしては
- グレースケール画像は0から1の範囲なので、255で割っている。
- 行列を可視化するような形に出力されるため、デフォルトではxが下方向、yが右方向になっている。これを変更するために、
imgを'で転置した。
です。Fortranと異なりJuliaはモダンな言語ですから、画像や動画の扱いが容易です。
万有引力に従う天体の運動
こちらのサイトでは、万有引力に従う粒子の運動の解説とFortranコードの説明があります。ここにあるFortranコードをJuliaに移植してみることにします。
基本方程式
二つの物体の間には物体の質量に比例しその距離の二乗に反比例する引力（万有引力）
\[F = - G \frac{Mm}{r^2}\]
が働きます。この力があるため、太陽の周りを地球などが回っているわけです。これをコードで示すことにします。
無次元化された運動方程式を
\[\dot{x} = u, \dot{u} = - \frac{x}{|{\bm r}|^3}\]
\[\dot{y} = v, \dot{v} = - \frac{y}{|{\bm r}|^3}\]
として、天体の動きを追いかけることとします。なお、簡単のためz方向には天体を固定し動かないこととしました。
方程式の時間発展にはいくつかの方法があります。微分方程式が
\[\dot{x} = u, \dot{u} = \Phi(x)\]
と書ける時の時間発展法について述べます。方程式のエネルギーを保存させるためにはシンプレクティック解法と呼ばれる手法を用いる必要がありまして、ここではLeap-frog法とVelocity Verlet法を紹介します。
Leap-frog法
\[u(t + \Delta t/2) = u(t - \Delta/2) + \Phi(x) \Delta t\]
\[x(t + \Delta t) = x(t) + u(t + \Delta t/2) \Delta t\]
という形になりまして、$\Delta/2$だけ速度uを動かしてからxを更新します。tでの速度は
\[u(t) = (u(t - \Delta/2) + u(t + \Delta/2))/2\]
と計算できます。最初のステップでは$\Delta(-\Delta t/2)$での値が必要ですので、これは方程式を逆に解いて
\[\Delta(-\Delta t/2) = u(0) - \Phi(x) \Delta /2\]
とします。
Velocity Verlet法
Velocity Verlet法の場合には、
\[x(t + \Delta t) = x(t) + u(t) \Delta t + \Phi(x) \frac{(\Delta t)^2}{2}\]
\[u(t + \Delta t) = u(t) + (\Phi(x(t)) + \Phi(x(t+\Delta t))) \frac{\Delta t}{2}\]
とします。
Leap-frog法による二体問題
Fortranで書かれたソースコードはこちらにあります。ここにあるコードをJuliaに移植してみました。
コードは
module Subfunc
export fx,fu,fy,fv
function fx(u,v)
return u
end
function fu(x,y)
r = sqrt(x^2+y^2)
xx = -x / r^3
return xx
end
function fy(u,v)
return v
end
function fv(x,y)
r = sqrt(x^2+y^2)
xx = -y / r^3
return xx
end
end
using .Subfunc
function celestial()
n=10^8
noutmax = 10000
x = 1.0
y = 0.0
u = 0.0
v = 1.0
t = 0.0
uold = 0.0
vold = 0.0
E0 = 0.0
u_b = 0.0
v_b = 0.0
tmax = 10.0e5
delta_t = (tmax-t)/n
fp = open("outputLF_julia.dat","w")
println(fp,"t_lf, x_lf, y_lf, E_lf, dE_lf")
if n > noutmax
dout = n ÷ noutmax
else
dout = 1
end
for i=1:n-1
if i == 1
E0 = 0.5*(u^2+v^2) - 1.0/sqrt(x^2+y^2)
u_b = u - 0.5*fu(x,y)*delta_t
v_b = v - 0.5*fv(x,y)*delta_t
uold = u
vold = v
end
u = uold
v = vold
E = 0.5*(u^2+v^2) - 1.0/sqrt(x^2+y^2)
dE = abs((E - E0)/E0)
if (i-1) % dout == 0
println(fp,"\t",t,"\t,",x,"\t,",y,"\t,",E,"\t,",dE)
end
u_a = u_b + fu(x,y)*delta_t
v_a = v_b + fv(x,y)*delta_t
x = x + u_a*delta_t
y = y + v_a*delta_t
uold = (u_a + u_b)/2.0
vold = (v_a + v_b)/2.0
u_b = u_a
v_b = v_a
t = t + delta_t
end
end
@time celestial()となります。リンク先のFortranコードを忠実に変換した形になっているかと思います。計算速度もほぼ同じになっているのが確認できると思います。
Velocity Verlet法によるN体問題
次に、N体問題を解いてみます。こちらのページを参考にしています。 複数の粒子がある場合は万有引力はそれぞれの粒子から受けますので、
\[F({\bm x}_i) = - G \sum_{j} \frac{m_i m_j}{|{\bm r}_i - {\bm r}_j|^2}\]
となります。
それでは、太陽、地球、金星、火星、そして月の5つの天体が存在する場合の軌道を調べてみましょう。そのFortran版のコードはこちらにあります。 このコードを手元のPCでコンパイルし、実行しますと、
current step: 100000 /total step: 1000000
current step: 200000 /total step: 1000000
current step: 300000 /total step: 1000000
current step: 400000 /total step: 1000000
current step: 500000 /total step: 1000000
current step: 600000 /total step: 1000000
current step: 700000 /total step: 1000000
current step: 800000 /total step: 1000000
current step: 900000 /total step: 1000000
Time step: 6.3717374890228697E-005
Average error of energy update: 1.0575727271732502E-021
CPU time: 0.38528799999999996 という結果が得られました。約0.4秒で実行が終わっていますね。このコードを素朴にJuliaのコードにしたものが
module Globals
export G,ndim,nbody
const G = 6.67384e-11
const ndim = 2
const nbody = 5
end
module Mod_Var_trajectory
using ..Globals
export Var_trajectory
Base.@kwdef struct Var_trajectory
coord::Array{Float64,2} = zeros(Float64,ndim,nbody)
veloc::Array{Float64,2} = zeros(Float64,ndim,nbody)
end
end
module Mod_Param_planet
using ..Globals
export Param_planet
Base.@kwdef struct Param_planet
mass::Array{Float64,1} = zeros(Float64,nbody)
end
end
module Cal_ene
using ..Globals
using ..Mod_Param_planet
export cal_kin,cal_pot
function cal_kin(param::Param_planet,v)
e = 0.0
for i=1:nbody
v2 = 0.0
for j=1:ndim
v2 += v[j,i]^2
end
e += 0.5*param.mass[i]*v2
end
return e
end
function cal_pot(param::Param_planet,c)
e = 0.0
rij = zeros(Float64,ndim)
for i=1:nbody-1
for j=i+1:nbody
rij[:] = c[:,j]-c[:,i]
rij2 = 0.0
for k=1:ndim
rij2 += rij[k]^2
end
rij_abs = sqrt(rij2)
e += param.mass[i]*param.mass[j]/rij_abs
end
end
return e
end
end
module Cal_force
using ..Globals
using ..Mod_Param_planet
export force!
function force!(param::Param_planet,c,a)
f = zeros(Float64,ndim,nbody)
rij = zeros(Float64,ndim)
fij = zeros(Float64,ndim)
for i = 1:nbody-1
for j = i+1:nbody
rij[:] = c[:,j]-c[:,i]
rij2 = 0.0
for k=1:ndim
rij2 += rij[k]^2
end
rij_abs = sqrt(rij2)
inv_rij_abs = 1.0/rij_abs
ftmp = param.mass[i]*param.mass[j]*(inv_rij_abs^3)
@. fij = ftmp*rij
f[:,i] = f[:,i] + fij[:]
f[:,j] = f[:,j] - fij[:]
end
end
for i = 1:nbody
a[:,i] = f[:,i] / param.mass[i]
end
end
end
using .Globals
using .Mod_Var_trajectory
using .Mod_Param_planet
using .Cal_ene
using .Cal_force
function celestial()
n=10^6
noutmax=10000
nlive=10
accel = zeros(Float64,ndim,nbody)
#!----- mass
param_planet = Param_planet()
mass = param_planet.mass
mass[1] = 1.989e30 # mass of sun[kg]
mass[2] = 5.972e24 # mass of earth[kg]
mass[3] = 4.869e24 # venus
mass[4] = 6.419e23 # mars
mass[5] = 7.348e22 # mass of moon[kg]
var_trajectory = Var_trajectory()
coord = var_trajectory.coord
#!----- initial conditions
coord[1,1] = 0.0
coord[2,1] = 0.0
coord[1,2] = 1.496e11 #!average radius of revolution of earth[m]
coord[2,2] = 0.0
coord[1,3] = 1.082e11
coord[2,3]= 0.0
coord[1,4] = 2.279e11
coord[2,4] = 0.0
coord[1,5] = 3.844e8 #!average radius of revolution of moon[m]
coord[1,5] += coord[1,2]
coord[2,5] = 0.0
veloc = var_trajectory.veloc
veloc[1,1] = 0.0
veloc[2,1] = 0.0
veloc[1,2] = 0.0
veloc[2,2] = 2.978e4 #!average velocity of earth[m/s]
veloc[1,3] = 0.0
veloc[2,3] = 3.502e4
veloc[1,4] = 0.0
veloc[2,4] = 2.413e4
veloc[1,5] = 0.0
veloc[2,5] = 1.022e3 #!average velocity of moon[m/s]
veloc[2,5] += veloc[2,2]
#!----- time
t = 0.0
tmax = 3.2e8 #!365day*24hour*60min*60sec = 31536000sec = 3.2d7
#!----- normalizse
cnor = sqrt((coord[1,2]-coord[1,1])^2 + (coord[2,2]-coord[2,1])^2)
mnor = mass[1]
tnor = sqrt(cnor^3/(G*mnor))
vnor = cnor/tnor
@. var_trajectory.coord = coord/cnor
@. var_trajectory.veloc = veloc/vnor
@. param_planet.mass = mass/mnor
t = t/tnor
tmax = tmax/tnor
dt = (tmax-t)/n
dt2 = 0.5*dt
sum_dE = 0.0
E0 = 0.0
Eold = 0.0
fp = open("outputNvV_julia_vv.dat","w")
println(fp,"t_NvV, E_NvV, dE_NvV, x1_NvV, y1_NvV, x2_NvV, y2_NvV,
x3_NvV, y3_NvV, x4_NvV, y4_NvV, x5_NvV, y5_NvV")
nliveout = n ÷ nlive
if n > noutmax
dout = n ÷ noutmax
else
dout = 1
end
for i = 1:n-1
if i % nliveout == 0
println("current step: $i,/total step: $n")
end
if i == 1
Ev = cal_kin(param_planet,veloc)
Ep = cal_pot(param_planet,coord)
E0 = Ev-Ep
Eold = E0
end
Ev = cal_kin(param_planet,veloc)
Ep = cal_pot(param_planet,coord)
E = Ev-Ep
dE = abs((E - E0)/E0)
sum_dE += abs(E-Eold)
if (i-1) % dout == 0
println(fp,
t,"\t",',',E,"\t",',',dE,"\t",',',coord[1,1],"\t",',',coord[2,1],"\t",',',coord[1,2],"\t",',',coord[2,2],"\t",',',coord[1,3]
,"\t",',',coord[2,3],"\t",',',coord[1,4],"\t",',',coord[2,4],"\t",',',coord[1,5],"\t",',',coord[2,5]
)
end
#println(coord)
if i == 1
force!(param_planet,coord,accel)
end
#println(accel)
#exit()
@. veloc += dt2*accel
@. coord += dt*veloc
force!(param_planet,coord,accel)
@. veloc += dt2*accel
t = t + dt
Eold = E
end
close(fp)
println("Time step: ",dt)
println("Average error of energy update: ",sum_dE/n)
end
@time celestial()こちらです。なるべく元のコードと同じ形になるように移植しました。 このコードのポイントは
- Fortranでのmoduleを用いた「グローバル変数」のような扱い（ここでは
mass等）を、全て独自型としてまとめて引数に入れるようにした - 1.を行った結果
var_trajectory.massのように変数名が長くなってしまったので、mass = var_trajectory.massのように参照用の変数を導入した。このようにするとmassの中身を書き換えたときにvar_trajectory.massも書き換わる。veloc = var_trajectory.velocやcoord = var_trajectory.coordも同様 - 配列の要素にまとめて入れるために
@.マクロを使用した。@. var_trajectory.coord = coord/cnorなどで使用
となります。このコードですが、同じPCで実行すると、
current step: 100000,/total step: 1000000
current step: 200000,/total step: 1000000
current step: 300000,/total step: 1000000
current step: 400000,/total step: 1000000
current step: 500000,/total step: 1000000
current step: 600000,/total step: 1000000
current step: 700000,/total step: 1000000
current step: 800000,/total step: 1000000
current step: 900000,/total step: 1000000
Time step: 6.37173748902287e-5
Average error of energy update: 1.0575727271732502e-21
5.327849 seconds (134.88 M allocations: 12.121 GiB, 4.46% gc time, 3.86% compilation time)となり、"Average error of energy update"は倍精度の範囲で結果が一致していることがわかります。ただし、計算時間が5.3秒もかかってしまっていて遅いです。
これを高速化する方法は簡単です。メモリーアロケーションが134.88 M allocationsとなっていますから、これは不必要なメモリ確保が沢山生じていることを意味しています。どこにそれがあるかをみてみますと、
rij[:] = c[:,j]-c[:,i]などです。右辺のc[:,j]は行列でいうところの1列だけを取り出しています。Juliaでは:を右辺で使って配列の一部範囲を取り出す場合、配列のコピーが発生してしまいます。これを避けるには@viewというマクロを使うか、素朴に
for k=1:length(rij)
rij[k] = c[k,j]-c[k,i]
endのようにループするかです。ループの方がシンプルなので、複数の該当箇所をループに置き換えてみると、
module Globals
export G,ndim,nbody
const G = 6.67384e-11
const ndim = 2
const nbody = 5
end
module Mod_Var_trajectory
using ..Globals
export Var_trajectory
Base.@kwdef struct Var_trajectory
coord::Array{Float64,2} = zeros(Float64,ndim,nbody)
veloc::Array{Float64,2} = zeros(Float64,ndim,nbody)
end
end
module Mod_Param_planet
using ..Globals
export Param_planet
Base.@kwdef struct Param_planet
mass::Array{Float64,1} = zeros(Float64,nbody)
end
end
module Cal_ene
using ..Globals
using ..Mod_Param_planet
export cal_kin,cal_pot
function cal_kin(param::Param_planet,v)
e = 0.0
for i=1:nbody
v2 = 0.0
for j=1:ndim
v2 += v[j,i]^2
end
e += 0.5*param.mass[i]*v2
end
return e
end
function cal_pot(param::Param_planet,c)
e = 0.0
rij = zeros(Float64,ndim)
for i=1:nbody-1
for j=i+1:nbody
for k=1:ndim
rij[k] = c[k,j]-c[k,i]
end
rij2 = 0.0
for k=1:ndim
rij2 += rij[k]^2
end
rij_abs = sqrt(rij2)
e += param.mass[i]*param.mass[j]/rij_abs
end
end
return e
end
end
module Cal_force
using ..Globals
using ..Mod_Param_planet
export force!
function force!(param::Param_planet,c,a)
f = zeros(Float64,ndim,nbody)
rij = zeros(Float64,ndim)
fij = zeros(Float64,ndim)
for i = 1:nbody-1
for j = i+1:nbody
for k=1:ndim
rij[k] = c[k,j]-c[k,i]
end
rij2 = 0.0
for k=1:ndim
rij2 += rij[k]^2
end
rij_abs = sqrt(rij2)
inv_rij_abs = 1.0/rij_abs
ftmp = param.mass[i]*param.mass[j]*(inv_rij_abs^3)
@. fij = ftmp*rij
for k=1:ndim
f[k,i] = f[k,i] + fij[k]
f[k,j] = f[k,j] - fij[k]
end
end
end
for i = 1:nbody
for k=1:ndim
a[k,i] = f[k,i] / param.mass[i]
end
end
end
end
using .Globals
using .Mod_Var_trajectory
using .Mod_Param_planet
using .Cal_ene
using .Cal_force
function celestial()
n=10^6
noutmax=10000
nlive=10
accel = zeros(Float64,ndim,nbody)
#!----- mass
param_planet = Param_planet()
mass = param_planet.mass
mass[1] = 1.989e30 # mass of sun[kg]
mass[2] = 5.972e24 # mass of earth[kg]
mass[3] = 4.869e24 # venus
mass[4] = 6.419e23 # mars
mass[5] = 7.348e22 # mass of moon[kg]
var_trajectory = Var_trajectory()
coord = var_trajectory.coord
#!----- initial conditions
coord[1,1] = 0.0
coord[2,1] = 0.0
coord[1,2] = 1.496e11 #!average radius of revolution of earth[m]
coord[2,2] = 0.0
coord[1,3] = 1.082e11
coord[2,3]= 0.0
coord[1,4] = 2.279e11
coord[2,4] = 0.0
coord[1,5] = 3.844e8 #!average radius of revolution of moon[m]
coord[1,5] += coord[1,2]
coord[2,5] = 0.0
veloc = var_trajectory.veloc
veloc[1,1] = 0.0
veloc[2,1] = 0.0
veloc[1,2] = 0.0
veloc[2,2] = 2.978e4 #!average velocity of earth[m/s]
veloc[1,3] = 0.0
veloc[2,3] = 3.502e4
veloc[1,4] = 0.0
veloc[2,4] = 2.413e4
veloc[1,5] = 0.0
veloc[2,5] = 1.022e3 #!average velocity of moon[m/s]
veloc[2,5] += veloc[2,2]
#!----- time
t = 0.0
tmax = 3.2e8 #!365day*24hour*60min*60sec = 31536000sec = 3.2d7
#!----- normalizse
cnor = sqrt((coord[1,2]-coord[1,1])^2 + (coord[2,2]-coord[2,1])^2)
mnor = mass[1]
tnor = sqrt(cnor^3/(G*mnor))
vnor = cnor/tnor
@. var_trajectory.coord = coord/cnor
@. var_trajectory.veloc = veloc/vnor
@. param_planet.mass = mass/mnor
t = t/tnor
tmax = tmax/tnor
dt = (tmax-t)/n
dt2 = 0.5*dt
sum_dE = 0.0
E0 = 0.0
Eold = 0.0
fp = open("outputNvV_julia_vv.dat","w")
println(fp,"t_NvV, E_NvV, dE_NvV, x1_NvV, y1_NvV, x2_NvV, y2_NvV,
x3_NvV, y3_NvV, x4_NvV, y4_NvV, x5_NvV, y5_NvV")
nliveout = n ÷ nlive
if n > noutmax
dout = n ÷ noutmax
else
dout = 1
end
for i = 1:n-1
if i % nliveout == 0
println("current step: $i,/total step: $n")
end
if i == 1
Ev = cal_kin(param_planet,veloc)
Ep = cal_pot(param_planet,coord)
E0 = Ev-Ep
Eold = E0
end
Ev = cal_kin(param_planet,veloc)
Ep = cal_pot(param_planet,coord)
E = Ev-Ep
dE = abs((E - E0)/E0)
sum_dE += abs(E-Eold)
if (i-1) % dout == 0
println(fp,
t,"\t",',',E,"\t",',',dE,"\t",',',coord[1,1],"\t",',',coord[2,1],"\t",',',coord[1,2],"\t",',',coord[2,2],"\t",',',coord[1,3]
,"\t",',',coord[2,3],"\t",',',coord[1,4],"\t",',',coord[2,4],"\t",',',coord[1,5],"\t",',',coord[2,5]
)
end
#println(coord)
if i == 1
force!(param_planet,coord,accel)
end
#println(accel)
#exit()
@. veloc += dt2*accel
@. coord += dt*veloc
force!(param_planet,coord,accel)
@. veloc += dt2*accel
t = t + dt
Eold = E
end
close(fp)
println("Time step: ",dt)
println("Average error of energy update: ",sum_dE/n)
end
@time celestial()となります。 この実行結果は、
current step: 100000,/total step: 1000000
current step: 200000,/total step: 1000000
current step: 300000,/total step: 1000000
current step: 400000,/total step: 1000000
current step: 500000,/total step: 1000000
current step: 600000,/total step: 1000000
current step: 700000,/total step: 1000000
current step: 800000,/total step: 1000000
current step: 900000,/total step: 1000000
Time step: 6.37173748902287e-5
Average error of energy update: 1.0575727271732502e-21
0.700156 seconds (4.88 M allocations: 509.987 MiB, 12.72% gc time, 30.04% compilation time)となり、計算速度がFortranとほぼ同じオーダーとなります。さらに、今は@timeで測定していますが、1回目はコードの最適化の時間が含まれてしまいます。ですので
@time celestial()
@time celestial()と二回実行してみますと、
current step: 100000,/total step: 1000000
current step: 200000,/total step: 1000000
current step: 300000,/total step: 1000000
current step: 400000,/total step: 1000000
current step: 500000,/total step: 1000000
current step: 600000,/total step: 1000000
current step: 700000,/total step: 1000000
current step: 800000,/total step: 1000000
current step: 900000,/total step: 1000000
Time step: 6.37173748902287e-5
Average error of energy update: 1.0575727271732502e-21
0.678261 seconds (4.88 M allocations: 509.987 MiB, 13.11% gc time, 30.31% compilation time)
current step: 100000,/total step: 1000000
current step: 200000,/total step: 1000000
current step: 300000,/total step: 1000000
current step: 400000,/total step: 1000000
current step: 500000,/total step: 1000000
current step: 600000,/total step: 1000000
current step: 700000,/total step: 1000000
current step: 800000,/total step: 1000000
current step: 900000,/total step: 1000000
Time step: 6.37173748902287e-5
Average error of energy update: 1.0575727271732502e-21
0.365792 seconds (4.42 M allocations: 486.000 MiB, 2.28% gc time)となり、二回目はFortranとほぼ同じスピードとなっていることがわかります。
Velocity Verlet法によるN体問題: Juliaっぽい書き方版
上で述べたコードはFortranコードを移植したものですから、すごく「Fortranっぽい」コードになっています。これを「Juliaっぽい」ものに書き換えてみようと思います。
まず、天体の情報を格納するstructを作りましょう。
struct Body
r::Array{Float64,1}
v::Array{Float64,1}
mass::Float64
f::Array{Float64,1}
end位置、速度、質量、そして力を格納します。 次に、天体たちの情報を格納するstructを作ります。
mutable struct Many_bodies
bodies::Array{Body,1}
ndim::Int64
nbody::Int64
rij_temp::Array{Float64,1}
fij_temp::Array{Float64,1}
function Many_bodies(ndim,nbody,mass)
bodies = Array{Body,1}(undef,nbody)
r = zeros(Float64,ndim)
for i=1:nbody
r = zero(r)
v = zero(r)
f = zero(r)
bodies[i] = Body(r,v,mass[i],f)
end
rij_temp = zero(r)
fij_temp = zero(r)
return new(bodies,ndim,nbody,rij_temp,fij_temp)
end
endここで、bodies::Array{Body,1}は、上で定義したBody型を要素に持つ配列です。Many_bodies(ndim,nbody,mass)はコンストラクタで、中身を初期化しています。つまり、引数が三つの関数Many_bodies(ndim,nbody,mass)を呼ぶことで、中身を指定しています。
次に、Many_bodies型からそれぞれの天体の情報を取り出すために
function Base.getindex(manybodies::Many_bodies,i)
return manybodies.bodies[i]
endを定義します。これは、manybodies[i]とmanybodiesを配列のように考えてその要素iを呼び出そうとするとどういう挙動をするのか、ということを決めています。引数がmanybodies::Many_bodiesとなっているのは、引数がMany_bodies型の時のみこの関数を呼ぶ、ということを指定しています。このように型ごとに関数を定義することで、オブジェクト指向プログラミングのような自由度を得ることができています。
さて、上のようなパーツを組み合わせてコードを作ります。その結果、
module Nbodys
export G,Many_bodies,Body
const G = 6.67384e-11
struct Body
r::Array{Float64,1}
v::Array{Float64,1}
mass::Float64
f::Array{Float64,1}
end
mutable struct Many_bodies
bodies::Array{Body,1}
ndim::Int64
nbody::Int64
rij_temp::Array{Float64,1}
fij_temp::Array{Float64,1}
function Many_bodies(ndim,nbody,mass)
bodies = Array{Body,1}(undef,nbody)
r = zeros(Float64,ndim)
for i=1:nbody
r = zero(r)
v = zero(r)
f = zero(r)
bodies[i] = Body(r,v,mass[i],f)
end
rij_temp = zero(r)
fij_temp = zero(r)
return new(bodies,ndim,nbody,rij_temp,fij_temp)
end
end
function Base.getindex(manybodies::Many_bodies,i)
return manybodies.bodies[i]
end
end
module Celestial_problem
using ..Nbodys
export cal_kin,cal_pot,cal_force!,update!
function cal_kin(manybodies::Many_bodies)
e = 0.0
for i=1:manybodies.nbody
v2 = 0.0
for j=1:manybodies.ndim
v2 += manybodies[i].v[j]^2
end
e += 0.5*manybodies[i].mass*v2
end
return e
end
function cal_pot(manybodies::Many_bodies)
ndim = manybodies.ndim
nbody = manybodies.nbody
e = 0.0
rij = manybodies.rij_temp
rij .= 0
for i=1:nbody-1
for j=i+1:nbody
for k=1:ndim
rij[k] = manybodies[j].r[k]-manybodies[i].r[k]
end
rij2 = 0.0
for k=1:ndim
rij2 += rij[k]^2
end
rij_abs = sqrt(rij2)
e += manybodies[i].mass*manybodies[j].mass/rij_abs
end
end
return e
end
function cal_force!(manybodies::Many_bodies)
ndim = manybodies.ndim
nbody = manybodies.nbody
rij = manybodies.rij_temp
rij .= 0
fij = manybodies.fij_temp
fij .= 0
for i=1:nbody
manybodies[i].f .= 0
end
for i = 1:nbody-1
for j = i+1:nbody
for k=1:ndim
rij[k] = manybodies[j].r[k]-manybodies[i].r[k]
end
rij2 = 0.0
for k=1:ndim
rij2 += rij[k]^2
end
rij_abs = sqrt(rij2)
inv_rij_abs = 1.0/rij_abs
ftmp = manybodies[i].mass*manybodies[j].mass*(inv_rij_abs^3)
@. fij = ftmp*rij
for k=1:ndim
manybodies[i].f[k] = manybodies[i].f[k] + fij[k]
manybodies[j].f[k] = manybodies[j].f[k] - fij[k]
end
end
end
end
function update!(manybodies::Many_bodies,dt,dt2)
nbody = manybodies.nbody
for k=1:nbody
@. manybodies[k].v += dt2*manybodies[k].f/manybodies[k].mass
@. manybodies[k].r += dt*manybodies[k].v
end
#exit()
cal_force!(manybodies)
for k=1:nbody
@. manybodies[k].v += dt2*manybodies[k].f/manybodies[k].mass
end
end
end
using .Celestial_problem
using .Nbodys
function celestial()
n=10^6
noutmax=10000
nlive=10
mass = zeros(Float64,5)
mass[1] = 1.989e30 # mass of sun[kg]
mass[2] = 5.972e24 # mass of earth[kg]
mass[3] = 4.869e24 # venus
mass[4] = 6.419e23 # mars
mass[5] = 7.348e22 # mass of moon[kg]]
mnor = mass[1]
@. mass = mass/mnor
ndim = 2
nbody = 5
manybodies = Many_bodies(ndim,nbody,mass)
manybodies[2].r[1] = 1.496e11 #!average radius of revolution of earth[m]
manybodies[3].r[1] = 1.082e11
manybodies[4].r[1] = 2.279e11
manybodies[5].r[1] = 3.844e8 #!average radius of revolution of moon[m]
manybodies[5].r[1] += manybodies[2].r[1]
manybodies[2].v[2] = 2.978e4 #!average velocity of earth[m/s]
manybodies[3].v[2] = 3.502e4
manybodies[4].v[2] = 2.413e4
manybodies[5].v[2] = 1.022e3 #!average velocity of moon[m/s]
manybodies[5].v[2] += manybodies[2].v[2]
#!----- time
t = 0.0
tmax = 3.2e8 #!365day*24hour*60min*60sec = 31536000sec = 3.2d7
#!----- normalizse
cnor = sqrt((manybodies[2].r[1]-manybodies[1].r[1])^2 +
(manybodies[2].r[2]-manybodies[1].r[2])^2)
tnor = sqrt(cnor^3/(G*mnor))
vnor = cnor/tnor
for i=1:nbody
@. manybodies[i].r = manybodies[i].r/cnor
@. manybodies[i].v = manybodies[i].v/vnor
end
t = t/tnor
tmax = tmax/tnor
dt = (tmax-t)/n
dt2 = 0.5*dt
sum_dE = 0.0
E0 = 0.0
Eold = 0.0
fp = open("outputNvV_julia_new.dat","w")
println(fp,"t_NvV, E_NvV, dE_NvV, x1_NvV, y1_NvV, x2_NvV, y2_NvV,
x3_NvV, y3_NvV, x4_NvV, y4_NvV, x5_NvV, y5_NvV")
nliveout = n ÷ nlive
if n > noutmax
dout = n ÷ noutmax
else
dout = 1
end
for i = 1:n-1
if i % nliveout == 0
println("current step: $i,/total step: $n")
end
if i == 1
Ev = cal_kin(manybodies)
Ep = cal_pot(manybodies)
E0 = Ev-Ep
Eold = E0
end
Ev = cal_kin(manybodies)
Ep = cal_pot(manybodies)
#exit()
E = Ev-Ep
dE = abs((E - E0)/E0)
sum_dE += abs(E-Eold)
if (i-1) % dout == 0
print(fp,t,"\t",',',E,"\t",',',dE,"\t",',')
for k=1:nbody
print(fp,manybodies[k].r[1],"\t",',',manybodies[k].r[2],"\t",',')
end
println(fp,"\t")
end
if i == 1
cal_force!(manybodies)
end
update!(manybodies,dt,dt2)
t = t + dt
Eold = E
end
close(fp)
println("Time step: ",dt)
println("Average error of energy update: ",sum_dE/n)
end
@time celestial()
@time celestial()
となりました。そして、実行すると、
current step: 100000,/total step: 1000000
current step: 200000,/total step: 1000000
current step: 300000,/total step: 1000000
current step: 400000,/total step: 1000000
current step: 500000,/total step: 1000000
current step: 600000,/total step: 1000000
current step: 700000,/total step: 1000000
current step: 800000,/total step: 1000000
current step: 900000,/total step: 1000000
Time step: 6.37173748902287e-5
Average error of energy update: 1.0582096959495853e-21
0.578360 seconds (848.32 k allocations: 79.538 MiB, 1.20% gc time, 35.65% compilation time)
current step: 100000,/total step: 1000000
current step: 200000,/total step: 1000000
current step: 300000,/total step: 1000000
current step: 400000,/total step: 1000000
current step: 500000,/total step: 1000000
current step: 600000,/total step: 1000000
current step: 700000,/total step: 1000000
current step: 800000,/total step: 1000000
current step: 900000,/total step: 1000000
Time step: 6.37173748902287e-5
Average error of energy update: 1.0582096959495853e-21
0.353715 seconds (390.21 k allocations: 55.551 MiB, 0.81% gc time)となり、以前のコードより若干良いパフォーマンスを持つコードとなりました。
Velocity Verlet法によるN体問題: さらに早く
上のコードは実はさらに高速化することができます。というのは、forループでndimやnbodyでループをしていますが、これは計算の最初に決まる数なので、定数になった方がコンパイラがうまく最適化してくれるからです。これを行うには、Many_bodiesをパラメトリック型Many_bodies{ndim,nbody}にします。パラメトリック型はこの型を宣言したときにndimとnbodyの値をパラメータとして取り込むことができます。これによって、Juliaはndimとnbodyをパラメータとして扱うことができ、高速化のためのループのアンローリングなどが可能となります。
mutable struct Many_bodies{ndim,nbody}
bodies::Array{Body,1}
rij_temp::Array{Float64,1}
fij_temp::Array{Float64,1}
function Many_bodies{ndim,nbody}(mass) where {ndim,nbody}
bodies = Array{Body,1}(undef,nbody)
r = zeros(Float64,ndim)
for i=1:nbody
r = zero(r)
v = zero(r)
f = zero(r)
bodies[i] = Body(r,v,mass[i],f)
end
rij_temp = zero(r)
fij_temp = zero(r)
return new{ndim,nbody}(bodies,rij_temp,fij_temp)
end
endこの型を使ったコードは以下の通りです。
module Nbodys
export G,Many_bodies,Body
const G = 6.67384e-11
struct Body
r::Array{Float64,1}
v::Array{Float64,1}
mass::Float64
f::Array{Float64,1}
end
mutable struct Many_bodies{ndim,nbody}
bodies::Array{Body,1}
rij_temp::Array{Float64,1}
fij_temp::Array{Float64,1}
function Many_bodies{ndim,nbody}(mass) where {ndim,nbody}
bodies = Array{Body,1}(undef,nbody)
r = zeros(Float64,ndim)
for i=1:nbody
r = zero(r)
v = zero(r)
f = zero(r)
bodies[i] = Body(r,v,mass[i],f)
end
rij_temp = zero(r)
fij_temp = zero(r)
return new{ndim,nbody}(bodies,rij_temp,fij_temp)
end
end
function Base.getindex(manybodies::Many_bodies,i)
return manybodies.bodies[i]
end
end
module Celestial_problem
using ..Nbodys
export cal_kin,cal_pot,cal_force!,update!
function cal_kin(manybodies::Many_bodies{ndim,nbody}) where {ndim,nbody}
e = 0.0
for i=1:nbody
v2 = 0.0
for j=1:ndim
v2 += manybodies[i].v[j]^2
end
e += 0.5*manybodies[i].mass*v2
end
return e
end
function cal_pot(manybodies::Many_bodies{ndim,nbody}) where {ndim,nbody}
e = 0.0
rij = manybodies.rij_temp
rij .= 0
for i=1:nbody-1
for j=i+1:nbody
for k=1:ndim
rij[k] = manybodies[j].r[k]-manybodies[i].r[k]
end
rij2 = 0.0
for k=1:ndim
rij2 += rij[k]^2
end
rij_abs = sqrt(rij2)
e += manybodies[i].mass*manybodies[j].mass/rij_abs
end
end
return e
end
function cal_force!(manybodies::Many_bodies{ndim,nbody}) where {ndim,nbody}
rij = manybodies.rij_temp
rij .= 0
fij = manybodies.fij_temp
fij .= 0
for i=1:nbody
manybodies[i].f .= 0
end
for i = 1:nbody-1
for j = i+1:nbody
for k=1:ndim
rij[k] = manybodies[j].r[k]-manybodies[i].r[k]
end
rij2 = 0.0
for k=1:ndim
rij2 += rij[k]^2
end
rij_abs = sqrt(rij2)
inv_rij_abs = 1.0/rij_abs
ftmp = manybodies[i].mass*manybodies[j].mass*(inv_rij_abs^3)
@. fij = ftmp*rij
for k=1:ndim
manybodies[i].f[k] = manybodies[i].f[k] + fij[k]
manybodies[j].f[k] = manybodies[j].f[k] - fij[k]
end
end
end
end
function update!(manybodies::Many_bodies{ndim,nbody},dt,dt2) where {ndim,nbody}
for k=1:nbody
@. manybodies[k].v += dt2*manybodies[k].f/manybodies[k].mass
@. manybodies[k].r += dt*manybodies[k].v
end
#exit()
cal_force!(manybodies)
for k=1:nbody
@. manybodies[k].v += dt2*manybodies[k].f/manybodies[k].mass
end
end
end
using .Celestial_problem
using .Nbodys
function ce2()
mass = zeros(Float64,5)
manybodies = Many_bodies{2,5}(mass)
end
function celestial()
n=10^6
noutmax=10000
nlive=10
mass = zeros(Float64,5)
mass[1] = 1.989e30 # mass of sun[kg]
mass[2] = 5.972e24 # mass of earth[kg]
mass[3] = 4.869e24 # venus
mass[4] = 6.419e23 # mars
mass[5] = 7.348e22 # mass of moon[kg]]
mnor = mass[1]
@. mass = mass/mnor
ndim = 2
nbody = 5
manybodies = Many_bodies{ndim,nbody}(mass)
manybodies[2].r[1] = 1.496e11 #!average radius of revolution of earth[m]
manybodies[3].r[1] = 1.082e11
manybodies[4].r[1] = 2.279e11
manybodies[5].r[1] = 3.844e8 #!average radius of revolution of moon[m]
manybodies[5].r[1] += manybodies[2].r[1]
manybodies[2].v[2] = 2.978e4 #!average velocity of earth[m/s]
manybodies[3].v[2] = 3.502e4
manybodies[4].v[2] = 2.413e4
manybodies[5].v[2] = 1.022e3 #!average velocity of moon[m/s]
manybodies[5].v[2] += manybodies[2].v[2]
#!----- time
t = 0.0
tmax = 3.2e8 #!365day*24hour*60min*60sec = 31536000sec = 3.2d7
#!----- normalizse
cnor = sqrt((manybodies[2].r[1]-manybodies[1].r[1])^2 +
(manybodies[2].r[2]-manybodies[1].r[2])^2)
tnor = sqrt(cnor^3/(G*mnor))
vnor = cnor/tnor
for i=1:nbody
@. manybodies[i].r = manybodies[i].r/cnor
@. manybodies[i].v = manybodies[i].v/vnor
end
t = t/tnor
tmax = tmax/tnor
dt = (tmax-t)/n
dt2 = 0.5*dt
sum_dE = 0.0
E0 = 0.0
Eold = 0.0
fp = open("outputNvV_julia_new_fastvinus.dat","w")
println(fp,"t_NvV, E_NvV, dE_NvV, x1_NvV, y1_NvV, x2_NvV, y2_NvV,
x3_NvV, y3_NvV, x4_NvV, y4_NvV, x5_NvV, y5_NvV")
nliveout = n ÷ nlive
if n > noutmax
dout = n ÷ noutmax
else
dout = 1
end
for i = 1:n-1
if i % nliveout == 0
println("current step: $i,/total step: $n")
end
if i == 1
Ev = cal_kin(manybodies)
Ep = cal_pot(manybodies)
E0 = Ev-Ep
Eold = E0
end
Ev = cal_kin(manybodies)
Ep = cal_pot(manybodies)
#exit()
E = Ev-Ep
dE = abs((E - E0)/E0)
sum_dE += abs(E-Eold)
if (i-1) % dout == 0
print(fp,t,"\t",',',E,"\t",',',dE,"\t",',')
for k=1:nbody
print(fp,manybodies[k].r[1],"\t",',',manybodies[k].r[2],"\t",',')
end
println(fp,"\t")
end
if i == 1
cal_force!(manybodies)
end
update!(manybodies,dt,dt2)
t = t + dt
Eold = E
end
close(fp)
println("Time step: ",dt)
println("Average error of energy update: ",sum_dE/n)
end
@time celestial()
@time celestial()
コードは型の宣言manybodies = Many_bodies{ndim,nbody}(mass)以外はndimとnbodyまわり以外全く変えていません。これだけの変更ですが、このコードを実行すると、
current step: 100000,/total step: 1000000
current step: 200000,/total step: 1000000
current step: 300000,/total step: 1000000
current step: 400000,/total step: 1000000
current step: 500000,/total step: 1000000
current step: 600000,/total step: 1000000
current step: 700000,/total step: 1000000
current step: 800000,/total step: 1000000
current step: 900000,/total step: 1000000
Time step: 6.37173748902287e-5
Average error of energy update: 1.0582096959495853e-21
0.540105 seconds (848.32 k allocations: 79.538 MiB, 1.89% gc time, 38.92% compilation time)
current step: 100000,/total step: 1000000
current step: 200000,/total step: 1000000
current step: 300000,/total step: 1000000
current step: 400000,/total step: 1000000
current step: 500000,/total step: 1000000
current step: 600000,/total step: 1000000
current step: 700000,/total step: 1000000
current step: 800000,/total step: 1000000
current step: 900000,/total step: 1000000
Time step: 6.37173748902287e-5
Average error of energy update: 1.0582096959495853e-21
0.315212 seconds (390.21 k allocations: 55.551 MiB, 0.99% gc time)となり、少し高速化します。これでgfortranの最適化オプションなしよりは速くなりました。
なお、gfortranのO3オプションをつけた場合の結果は
current step: 100000 /total step: 1000000
current step: 200000 /total step: 1000000
current step: 300000 /total step: 1000000
current step: 400000 /total step: 1000000
current step: 500000 /total step: 1000000
current step: 600000 /total step: 1000000
current step: 700000 /total step: 1000000
current step: 800000 /total step: 1000000
current step: 900000 /total step: 1000000
Time step: 6.3717374890228697E-005
Average error of energy update: 1.0521106352128808E-021
CPU time: 0.16022700000000001 ですので、まだ負けています。現在、Juliaのコードをさらに高速化する方法について試行錯誤しています。
2次元Ising模型
次は、2次元Ising模型のFORTRAN77コードをJuliaに移植してみましょう。
理論的背景
2次元Ising模型
古典スピン系であるIsing(イジング)模型のハミルトニアンは
\[H = -J \sum_{\langle i j \rangle} \sigma_i \sigma_j - h \sum_i \sigma_i\]
です。第二項は磁場の効果です。ここで、$\langle i j \rangle$は、最隣接格子点のみで和を取ることを意味していて、一次元系であれば、$j=i+1$などです。$\sigma_i$は$i$番目の格子点のスピンを表し、$+1$か$-1$を取ります。 統計力学において、物理量$A$の期待値は
\[\langle A \rangle = \frac{1}{Z} \sum_{\cal C} \left[ \exp \left( -\frac{H({\cal C})}{k_{\rm B} T} \right) A({\cal C}) \right]\]
と書けます。ここで、$H({\cal C})$は、あるスピン配置${\cal C}$でのハミルトニアン、$A({\cal C})$はその時の物理量$A$の値です。 $k_{\rm B}$はボルツマン定数、$T$は温度です。 $Z$は分配関数であり、
\[Z = \sum_{\cal C}\exp \left( -\frac{H({\cal C})}{k_{\rm B} T} \right)\]
と定義されています。 つまり、すべての可能なスピン配置に関して和を取れば、物理量が計算できます。 すべての可能なスピン配置の数${\cal N}$は、$N$個の格子点を持つ系の場合、各サイトで$-1$か$1$の二通りの状態を取れるので、
\[{\cal N} = 2^N\]
です。
\[L_x \times L_y\]
の正方格子の2次元Ising模型を考えてみましょう。この時、あるサイトを${\bf i} = (i_x,i_y)$とすると、その最近接格子は、
\[{\bf d}_1 = (i_x+1,i_y),{\bf d}_2 = (i_x-1,i_y),{\bf d}_3 = (i_x,i_y+1), {\bf d}_4 = (i_x,i_y-1)\]
の4点です。この時、Ising模型は
\[H = -\frac{J}{2} \sum_{\bf i}^{L_x L_y} \sum_{l=1}^4 \sigma_{\bf i} \sigma_{{\bf i}}{}_{ + {\bf d}_l}- h \sum_{\bf i} \sigma_{\bf i}\]
となります。ここで、本来一つしかない$\sigma_1 \sigma_2$を$\sigma_1 \sigma_2 = (\sigma_1 \sigma_2 + \sigma_2 \sigma_1)/2$と分離したため、$1/2$の因子が先頭につきました。 そしてこれは、
\[H = -\frac{J}{2} \sum_{\bf i}^{L_x L_y} \sigma_{\bf i} \sum_{l=1}^4 \sigma_{{\bf i}}{}_{ + {\bf d}_l}- h \sum_{\bf i} \sigma_{\bf i}\]
\[= -\frac{J}{2} \sum_{\bf i}^{L_x L_y} \sigma_{\bf i} S_i - h \sum_{\bf i} \sigma_{\bf i}\]
\[S_i = \sum_{l=1}^4 \sigma_{{\bf i}}{}_{ + {\bf d}_l}\]
と書くことができます。よって、あるサイト${\bf i}$の隣接格子点におけるスピンの和$S_i$がそれぞれわかれば、全エネルギー$H$を計算できます。 あるサイト${\bf i}$の一つのスピンをフリップ$(\sigma_{\bf i} \rightarrow -\sigma_{\bf i})$した時、そのエネルギー差$\Delta E$は
\[\Delta E = 2J \sigma_{\bf i} S_i+ 2h \sigma_{\bf i}\]
となります。
メトロポリス法
あるサイトをランダムに選び、そのスピンをフリップさせることでスピン配置を更新するとします。この時、全サイト数を$N = L_x \times L_y$とすると、確率$1/N$でサイト数を選ぶ。そして、このようなスピンフリップであればプロセスは対称であるので、メトロポリス法でのある配置$C_i$から$C_j$への採択率は
\[A(C_i \rightarrow C_j) = {\rm min} \left(1, \frac{P(C_j)}{P(C_i)}\right)\]
となります。ここで、$P(C_i)$をボルツマン重み
\[P(C_i) = \exp \left( -\frac{H({\cal C}_i)}{k_{\rm B} T} \right)\]
と選べば、物理量を重みつきモンテカルロ法で計算できます。 スピン配置${\cal C}_k$のあるサイト${\bf i}$のスピン$\sigma_{\bf i}$をフリップさせスピン配置${\cal C}_k'$とする時、採択率に現れる重みの比は
\[\frac{P(C_k')}{P(C_k)} = \exp \left( -\frac{(H({\cal C}_k)-H({\cal C}_k'))}{k_{\rm B} T} \right)\]
\[= \exp \left( -\frac{\Delta E({\cal C}_k,{\bf i})}{k_{\rm B} T}\right)\]
となります。
FORTRAN77コード
G. パリージ、場の理論―統計論的アプローチ（吉岡書店）1993の20章には、2次元Ising模型のモンテカルロシミュレーションコードが記載されています。原著の出版年は1988年ですから、このコードはFortran90ではなく、FORTRAN77で書かれています。 このコードを丸写ししても現代のFortranコンパイラではうまく動きません(RAN関数の動作がおかしい)。ですので、2021年のgfortranでも動くようにコードを少し手直ししました。手直しして動くようになったコードがこちらです(微修正に過ぎませんので、教科書の元のコードの引用になるかと思います)。
PARAMETER maximal_side=100
INTEGER forward(maximal_side),backward(maximal_side)
COMMON /boundaries/ forward,backward
INTEGER side,spin(maximal_side,maximal_side)
INTEGER random_seed
REAL magnetic_field,magnetization_density
CALL read_input(side,number_of_iterations,beta,magnetic_field)
CALL get_random_seed(random_seed)
call srand(random_seed)
CALL compute_backward_and_forward(side)
CALL set_spin_to_1(spin,side)
DO 1 iteration=1,number_of_iterations
CALL one_Monte_Carlo_cycle(spin,side,random_seed,
1 beta,magnetic_field,energy_density,magnetization_density)
CALL write_output(iteration,energy_density,magnetization_
1 density)
1 CONTINUE
END
SUBROUTINE compute_backward_and_forward(side)
PARAMETER maximal_side=100
INTEGER forward(maximal_side),backward(maximal_side)
COMMON /boundaries/ forward,backward
INTEGER position,side
DO 1 position=1,side
forward(position)=MOD(position,side)+1
backward(position)=MOD((position-2+side),side)+1
1 CONTINUE
RETURN
END
SUBROUTINE set_spin_to_1(spin,side)
PARAMETER maximal_side=100
INTEGER side,spin(maximal_side,maximal_side)
INTEGER x_position,y_position
DO 1 x_position=1,side
DO 2 y_position=1,side
spin(x_position,y_position)=1
2 CONTINUE
1 CONTINUE
RETURN
END
SUBROUTINE one_Monte_Carlo_cycle(spin,side,random_seed,
1beta,magnetic_field,energy_density,magnetization_density)
PARAMETER maximal_side=100
INTEGER forward(maximal_side),backward(maximal_side)
COMMON /boundaries/ forward,backward
INTEGER side,random_seed,spin(maximal_side,maximal_side)
REAL magnetic_field,magnetization_density
INTEGER current_spin,sum_of_the_neighbours
INTEGER x_position,y_position
total_magnetization=0
total_energy=0
DO 5 x_position=1,side
DO 6 y_position=1,side
current_spin=spin(x_position,y_position)
sum_of_the_neighbours=
1 spin(x_position,forward(y_position))+
2 spin(x_position,backward(y_position))+
3 spin(forward(x_position),y_position)+
4 spin(backward(x_position),y_position)
effective_force=sum_of_the_neighbours+magnetic_field
IF(EXP(-beta*effective_force*current_spin*2.)
1 .GT. RAND() ) THEN
new_spin = -current_spin
ELSE
new_spin=current_spin
END IF
spin(x_position,y_position)=new_spin
total_magnetization=total_magnetization+new_spin
total_energy = total_energy+
1 (.5*sum_of_the_neighbours+magnetic_field)*new_spin
6 CONTINUE
5 CONTINUE
magnetization_density=total_magnetization/FLOAT(side**2)
energy_density=total_energy/FLOAT(side**2)
RETURN
END
SUBROUTINE read_input(side,number_of_iterations,beta,
1magnetic_field)
REAL magnetic_field
INTEGER side
PARAMETER maximal_side=100
1 WRITE(6,*) 'Which is the length of the side ?
1 (Please less than',maximal_side,')'
READ(5,*) side
IF(side .GT. maximal_side .OR. side .LT. 1) THEN
WRITE(6,*) 'The value of side is not good'
go to 1
ENDIF
WRITE(6,*) 'How many iterations?'
READ(5,*) number_of_iterations
WRITE(6,*) 'Which is the value of beta?'
READ(5,*) beta
WRITE(6,*) 'Which is the value of the magnetic field?'
READ(5,*) magnetic_field
RETURN
END
SUBROUTINE get_random_seed(random_seed)
INTEGER random_seed
WRITE(6,*) 'Which is the random seed?'
WRITE(6,*) 'Please insert a positive odd number of 7-8 digits'
READ(5,*) random_seed
RETURN
END
SUBROUTINE write_output(iteration,energy_density,
1 magnetization_density)
REAL magnetization_density
WRITE(6,*) 'Iteration= ',iteration
WRITE(6,*) 'Energy density = ',energy_density
WRITE(6,*) 'Magnetization density= ',magnetization_density
RETURN
END
goto文やcontinue文、common文が使われています。program mainというFortran90でお馴染みのものもありません。そして、実数が全て単精度になっています。betaが温度の逆数になっていまして、$\beta_c = \tanh^{-1}(\sqrt(2) -1 ) = 0.44068679350977163$が相転移温度です。この$\beta_c$よりbetaが小さければ高温側、大きければ低温側となり、1サイトあたりの磁化であるMagnetization densityの値は高温で0、低温1となります。
このコードをJuliaに移植してみます。以下がそのコードです。
struct Boundaries
forward::Array{Int64,1}
backward::Array{Int64,1}
function Boundaries(side)
return new(zeros(Int64,side),zeros(Int64,side))
end
end
function compute_backward_and_forward(side)
boundaries = Boundaries(side)
forward = boundaries.forward
backward = boundaries.backward
for position=1:side
forward[position]=position % side +1
backward[position]=(position-2+side) % side+1
end
return boundaries
end
function set_spin_to_1!(spin,side)
spin[1:side,1:side] .= 1
end
function one_Monte_Carlo_cycle!(boundaries,spin,side,
beta,magnetic_field)
forward = boundaries.forward
backward = boundaries.backward
total_magnetization=0.0
total_energy=0.0
for x_position=1:side
for y_position=1:side
current_spin=spin[x_position,y_position]
sum_of_the_neighbours=
spin[x_position,forward[y_position]] +
spin[x_position,backward[y_position]]+
spin[forward[x_position],y_position]+
spin[backward[x_position],y_position]
effective_force=sum_of_the_neighbours+magnetic_field
if exp(-beta*effective_force*current_spin*2) > rand()
new_spin = -current_spin
else
new_spin = current_spin
end
spin[x_position,y_position]=new_spin
total_magnetization=total_magnetization+new_spin
total_energy = total_energy+
(0.5*sum_of_the_neighbours+magnetic_field)*new_spin
end
end
magnetization_density=total_magnetization/(side^2)
energy_density=total_energy/(side^2)
return energy_density,magnetization_density
end
function read_input(;maximal_side=100)
side = 0
while true
println("Which is the length of the side ? (Please less than ",maximal_side,")")
side = parse(Int64,readline())
if side > maximal_side || side < 1
println("The value of side is not good")
else
break
end
end
println("How many iterations?")
number_of_iterations = parse(Int64,readline())
println("Which is the value of beta?")
beta = parse(Float64,readline())
println("Which is the value of the magnetic field?")
magnetic_field = parse(Float64,readline())
return side,number_of_iterations,beta,magnetic_field
end
function get_random_seed()
println("Which is the random seed?")
println("Please insert a positive odd number of 7-8 digits")
random_seed = parse(Int64,readline())
return random_seed
end
function write_output(iteration,energy_density,
magnetization_density)
println("Iteration= ",iteration)
println("Energy density = ",energy_density)
println("Magnetization density= ",magnetization_density)
return
end
using Random
function main()
side,number_of_iterations,beta,magnetic_field = read_input(maximal_side=100)
spin = zeros(Float64,side,side)
random_seed = get_random_seed()
Random.seed!(random_seed)
boundaries = compute_backward_and_forward(side)
set_spin_to_1!(spin,side)
for iteration=1:number_of_iterations
energy_density,magnetization_density =
one_Monte_Carlo_cycle!(boundaries,spin,side,beta,magnetic_field)
write_output(iteration,energy_density,magnetization_density)
end
end
main()となります。
ポイントは、
- common文はstructで扱うようにした。これで不必要に引数が増えることがなくなり、どの関数がどの値を扱っているかがわかりやすくなった
- 現代では単精度でCPU上で計算する意味はもはやないので、全て倍精度実数とした
- goto文での無限ループの部分は
while trueで回すようにした - キーボードからの入力には
readline()を使った
です。
FortranコードをJuliaに移植するためのチェックシート
「Leap-frog法による二体問題」のFortranコードをJuliaに移植したときに意識した点をチェックシートにしました。何か他に気がついた点があれば適宜足していこうとは思っています。
- [ ] Fortranでの変数の定義は基本使わないので消す
- [ ]
end doのdoなどの削除：end doやend if、end functionなどは全部endに - [ ] doループ：
do i=1,3のようなものはfor i=1:3に - [ ] programの置き換え：
program mainのようなものはfunction main()として定義の次の行にmain()を書く - [ ] べき乗の書き換え：
**を一括置換で^に - [ ] 余り計算の書き換え：
mod(a,b)はa % bに - [ ]
dble(n)のdbleの除去： Juliaではnで問題なし - [ ] If文の書き換え：
if (a .eq. b) thenはif a == bに - [ ] 数字の書き換え： 倍精度の表現
d0などは一括置換でスペースに - [ ] 桁数の表現の書き換え： 倍精度の
d5などはe5とdをeに - [ ] forループ内の変数に注意：Juliaだとforループ内でのみ定義されている変数は外に出せないので、forループの外で
a = 0.0のように定義しておく - [ ] ファイルオープン関連：
open(10,filename="hoge.dat")などはfp = open("hoge.dat","w")に - [ ] 出力関連：
write(11,*) a,b,cなどはprintln(fp,"$a \t $b \t $c")に - [ ] module関連：FortranのmoduleはJuliaでもそのままmoduleに。ただし名前の最初を大文字にした方がJuliaの慣例に合う
- [ ] moduleの使用関連： Fortranで関数内で
use hogeとしている場合、Juliaではその関数の外でusing .hogeとして呼び出す - [ ] module内関数の使用： Fortranで
use hogeで使っている関数fは、Juliaではmodule内でexport fをした状態でusing .hogeを行う
こちらに加えて、Fortranが配列を使っている場合は以下のことも気にした方が良いと思います。
- [ ] 配列の定義の書き換え：
real(8) a(2,3)のようなものはa = zeros(Float64,2,3)などと書き換え。integer b(4,5)のような整数の場合にはb = zeros(Int64,4,5)とする - [ ]
allocatable属性の書き換え： Fortranでallocateされている場所でJuliaでは上のように配列を定義しておく - [ ] 配列の引数範囲が
real(8) a(-1:3,-3:4)のように1始まりではない場合：Juliaではusing OffsetArraysをしてから、a = OffsetArray(zeros(Float64,5,8), -1:3, -3:4)のようにする
グローバル変数関連は、以下のようにすると良いでしょう。
- [ ] common文はstructにしてしまう：common文は同じ並びのstructを作ってしまって、そのstructを関数の第一引数に入れるような形で取り回すと変数の変更の有無がわかりやすい
- [ ] module文内にあるsave属性の変数（グローバル変数のような変数）：common文の処方箋と同じ
- [ ] module文内にある変更されない変数：module内で
const a = 3のようにconstをつける