神戸大学 大学院システム情報学研究科 計算科学専攻 陰山 聡
【目次】
pause -1と書く必要があった。
前回の宿題は、「thermal_diffusion.f90を理解しておくこと」であった。 ここではポイントだけ説明する。
| モジュール名 | 内容 |
| constants | 定数 |
| parallel | 並列化 |
| temperature | 温度場 |
データ出力ルーチンを除くとコードの構造は以下のとおり：
call parallel__initialize
call temperature__initialize
do loop = 1 , LOOP_MAX
call temperature__update
end do
call temperature__finalize
call parallel__finalize
ルーチン名を見れば意味は明白であろう。
call parallel__initialize call temperature__initialize call temperature__output_1d_profile call temperature__output_2d_profile
に見られるように、関数名（サブルーチン名）は、「モジュール名」＋「アンダースコア二つ」＋「その関数の中身を示す言葉」で統一されている。 その関数がどのモジュールで定義されているかがすぐに分かるようにするために工夫である。
今考えている正方形領域は 一様かつ一定に加熱されている。 その加熱率 heat_source を（ある規格化の下で）4とすると、 定常状態における温度T(x,y)は以下のポアッソン方程式を満たす。
この式を丁寧に書くと、
である。
連続場としての温度場を x、y方向に一様な格子間隔（2次元カーテシアン格子）で離散化し、 2次精度の中心差分法で上のラプラシアンを近似すると、
(T(i+1,j)-2*T(i,j)+T(i-1,j))/dx^2 + (T(i,j+1)-2*T(i,j)+T(i,j-1))/dy^2 + heat_source = 0.
となる。ここでdxとdyはそれぞれx方向とy方向の格子間隔である。 特にdx=dyならば、その値をhと書いて、変形すると
T(i,j) = (T(i+1,j)+T(i-1,j)+T(i,j+1)+T(i,j-1))/4+heat_source*h^2/4.
が得られる。 ヤコビ法ではこの式が成り立つまで温度場Tを繰り返し計算する。
4*h^2/4
となる。この項は一定なので、ヤコビ法の繰り返し計算で毎回計算するのは馬鹿らしい。 そこで、thermal_diffusion.f90コードでは、subroutine temperature__initialize において
subroutine temperature__initialize real(DP) :: heat_source = 4.0 . . . heat = (heat_source/4) * h * h
としてあらかじめ（一度だけ）計算しておき、その値をheatという変数に保存した上で、 ヤコビ法の繰り返し部分（subroutine temperature__update）では、
subroutine temperature__update
integer :: i, j
.
.
.
do j = jstart, jend
do i = 1, MESH_SIZE
! work(i,j) = (temp(i-1,j)+temp(i+1,j)+temp(i,j-1)+temp(i,j+1))/4.0_DP &
! + heat_source*h*h ! Operation time is saved in below.
work(i,j) = (temp(i-1,j)+temp(i+1,j)+temp(i,j-1)+temp(i,j+1))*0.25_DP &
+ heat
end do
end do
のようにheatという定数を足すだけで済むようにしている。
今回は2次元データの可視化を行う。 前回と同様gnuplotを利用する。
gnuplotで2次元データをプロットする場合、そのデータのフォーマットは
x00 y00 関数値 x01 y00 関数値 x02 y00 関数値 ・ ・ ・ x09 y00 関数値 （空行） x00 y01 関数値 x01 y01 関数値 x02 y01 関数値 ・ ・ ・ x09 y01 関数値 （空行） ・ ・ ・ x09 y09 関数値
である。
thermal_diffusion.f90による2次元データ出力ルーチン temperature__output_2d_profile は、正方形上（x,y平面上）に分布する温度をすべて書き出すものである。 中身を見てみよう。
subroutine temperature__output_2d_profile
real(DP), dimension(0:MESH_SIZE+1, &
0:MESH_SIZE+1) :: prof
integer :: counter = 0 ! saved
integer :: ierr ! use for MPI
integer :: jstart_, jend_
character(len=4) :: serial_num ! put on file name
character(len=*), parameter :: base = "../data/temp.2d."
integer :: i, j
call set_jstart_and_jend(jstart_,jend_)
write(serial_num,'(i4.4)') counter
prof(:,:) = parallel__set_prof_2d(jstart_,jend_,temp(:,jstart_:jend_))
if ( myrank==0 ) then
open(10,file=base//serial_num)
do j = 0 , MESH_SIZE+1
do i = 0 , MESH_SIZE+1
write(10,*) i, j, prof(i,j)
end do
write(10,*)' ' ! gnuplot requires a blank line here.
end do
close(10)
end if
counter = counter + 1
end subroutine temperature__output_2d_profile
dataディレクトリ中の連番つきファイルtemp.2d.????の中身は以下のようになっているはず。 確認せよ。
55 35 0.1165598588705999
56 35 9.9624877672293416E-002
57 35 8.1734108631726782E-002
58 35 6.2857224006520482E-002
59 35 4.2963409431420671E-002
60 35 2.2021479795155254E-002
61 35 0.000000000000000
0 36 0.000000000000000
1 36 2.1867122785152873E-002
2 36 4.2655284590767971E-002
3 36 6.2396502500601705E-002
4 36 8.1122529495226178E-002
・
・
・
57 36 8.1122529495226178E-002
58 36 6.2396502500601705E-002
59 36 4.2655284590767971E-002
60 36 2.1867122785152873E-002
61 36 0.000000000000000
0 37 0.000000000000000
1 37 2.1680615643577435E-002
2 37 4.2282992534786040E-002
3 37 6.1839860798781864E-002
4 37 8.0383668447763873E-002
5 37 9.7946446090882752E-002
6 37 0.1145596782417825
・
・
・
# # a sample gnuplot script: sample_contour_lnes.gp # # [ line contours ] # set size square # same side lengths for x and y set xlabel 'i' # x-axis set ylabel 'j' # y-axis set xrange[0:62] # i-grid min & max set yrange[0:62] # j-grid min & max set nosurface # do not show surface plot unset ztics # do not show z-tics set contour base # enables contour lines set cntrparam levels 10 # draw 10 contours set view 0,0 # view from the due north set title 'Temperature at 100' splot '../data/temp.2d.0100' using 1:2:3 w l # with lines pause -1
!
! anime_contour_lines_generator.f90
!
! Usage:
! (1) check the values of mesh size m and counter_end.
! (2) pg95 this_code
! (3) ./a.out > anyname
! (4) gnuplot anyname
!
program anime_contour_lines_generator
implicit none
integer, parameter :: m = 61 ! mesh size
integer, parameter :: counter_end = 100
integer :: counter
character(len=*), parameter :: base='../data/temp.2d.'
character(len=4) :: serial_num
print *, "#"
print *, "# gnuplot script generated by anime_contour_lines_generator.f90"
print *, "#"
print *, "set size square # same side lengths for x and y"
print *, "set xlabel 'i' # x-axis"
print *, "set ylabel 'j' # y-axis"
print *, "set xrange[0:", m+1, "] # i-grid min & max"
print *, "set yrange[0:", m+1, "] # j-grid min & max"
print *, "set nosurface # do not show surface plot"
print *, "unset ztics # do not show z-tics"
print *, "set contour base # enables contour lines"
print *, "set cntrparam levels 10 # draw 10 contours"
print *, "set view 0,0 # view from the due north"
do counter = 0 , counter_end
write(serial_num,'(i4.4)') counter
print *, "set title 'Temperature at ", counter, "'"
print *, "splot '"//base//serial_num//"' using 1:2:3 w l"
select case (counter)
case (0)
print *, "pause 5"
case (counter_end)
print *, "pause -1"
case default
print *, "pause 0.2"
end select
end do
end program anime_contour_lines_generator
# # a sample gnuplot script: sample_contour_colors.gp # # [ color contours ] # set size square # same side lengths for x and y set xlabel 'i' # x-axis set ylabel 'j' # y-axis set xrange[0:62] # i-grid min & max set yrange[0:62] # j-grid min & max set palette defined (0 'blue', 0.15 'red', 0.3 'yellow') set nosurface # do not show surface plot unset ztics # do not show z-tics set pm3d at b # draw with colored contour set view 0,0 # view from the due north set title 'Temperature at 100 ' splot '../data/temp.2d.0100' using 1:2:3 pause -1
!
! anime_contour_colors_generator.f90
!
! Usage:
! (1) check the values of mesh size m and counter_end.
! (2) pg95 this_code
! (3) ./a.out > anyname
! (4) gnuplot anyname
!
program anime_contour_colors_generator
implicit none
integer, parameter :: m = 61 ! mesh size
integer, parameter :: counter_end = 100
integer :: counter
character(len=*), parameter :: base='../data/temp.2d.'
character(len=4) :: serial_num
print *, "#"
print *, "# gnuplot script generated by anime_contour_colors_generator.f90"
print *, "#"
print *, "set size square # same side lengths for x and y"
print *, "set xlabel 'i' # x-axis"
print *, "set ylabel 'j' # y-axis"
print *, "set xrange[0:", m+1, "] # i-grid min & max"
print *, "set yrange[0:", m+1, "] # j-grid min & max"
print *, "set palette defined (0 'blue', 0.15 'red', 0.3 'yellow')"
print *, "set nosurface # do not show surface plot"
print *, "unset ztics # do not show z-tics"
print *, "set pm3d at b # draw with colored contour"
print *, "set view 0,0 # view from the due north"
do counter = 0 , counter_end
write(serial_num,'(i4.4)') counter
print *, "set title 'Temperature at ", counter, "'"
print *, "splot '"//base//serial_num//"' using 1:2:3 "
select case (counter)
case (0)
print *, "pause 5"
case (counter_end)
print *, "pause -1"
case default
print *, "pause 0.2"
end select
end do
end program anime_contour_colors_generator
# # a sample gnuplot script: sample_birdseyeview.gp # # [ Bird's Eye View ] # set size square # same side lengths for x and y set xlabel 'i' # x-axis set ylabel 'j' # y-axis set xrange[0:62] # i-grid min & max set yrange[0:62] # j-grid min & max set contour base # enables contour lines set cntrparam levels 10 # draw 10 contours # set palette defined (0 'blue', 0.15 'red', 0.3 'yellow') # set pm3d # draw with colored contour set title 'Temperature at 100 ' splot '../data/temp.2d.0100' using 1:2:3 w l pause -1
!
! anime_birdseyeview_generator.f90
!
! Usage:
! (1) check the values of mesh size m and counter_end.
! (2) pg95 this_code
! (3) ./a.out > anyname
! (4) gnuplot anyname
!
program anime_birdseyeview_generator
implicit none
integer, parameter :: m = 61 ! mesh size
integer, parameter :: counter_end = 100
integer :: counter
character(len=*), parameter :: base='../data/temp.2d.'
character(len=4) :: serial_num
print *, "#"
print *, "# gnuplot script generated by anime_birdseyeview_generator.f90"
print *, "#"
print *, "set size square # same side lengths for x and y"
print *, "set xlabel 'i' # x-axis"
print *, "set ylabel 'j' # y-axis"
print *, "set xrange[0:", m+1, "] # i-grid min & max"
print *, "set yrange[0:", m+1, "] # j-grid min & max"
print *, "set zrange[-0.35:0.35] # temperature min & max"
print *, "set contour base # enables contour lines"
print *, "set cntrparam levels 10 # draw 10 contours"
do counter = 0 , counter_end
write(serial_num,'(i4.4)') counter
print *, "set title 'Temperature at ", counter, "'"
print *, "splot '"//base//serial_num//"' using 1:2:3 w l"
select case (counter)
case (0)
print *, "pause 5"
case (counter_end)
print *, "pause -1"
case default
print *, "pause 0.2"
end select
end do
end program anime_birdseyeview_generator
!
! anime_birdseyeview_rotation_generator.f90
!
! Usage:
! (1) check the values of mesh size m and counter_end.
! (2) pg95 this_code
! (3) ./a.out > anyname
! (4) gnuplot anyname
!
program anime_birdseyeview_rotation_generator
implicit none
integer, parameter :: m = 61 ! mesh size
integer, parameter :: counter_target = 100
integer, parameter :: angle_start = 30
integer, parameter :: angle_end = angle_start + 180
integer :: angle
character(len=*), parameter :: base='../data/temp.2d.'
character(len=4) :: serial_num
print *, "#"
print *, "# gnuplot script generated by anime_birdseyeview_generator.f90"
print *, "#"
print *, "set size square # same side lengths for x and y"
print *, "set xlabel 'i' # x-axis"
print *, "set ylabel 'j' # y-axis"
print *, "set xrange[1:", m+1, "] # i-grid min & max"
print *, "set yrange[1:", m+1, "] # j-grid min & max"
print *, "set zrange[-0.35:0.35] # temperature min & max"
print *, "set contour base # enables contour lines"
print *, "set cntrparam levels 10 # draw 10 contours"
do angle = angle_start , angle_end
write(serial_num,'(i4.4)') counter_target
print *, "set title 'Temperature at ", counter_target, "'"
print *, "set view 60, ", angle
print *, "splot '"//base//serial_num//"' using 1:2:3 w l"
select case (angle)
case (angle_start)
print *, "pause 5"
case default
print *, "pause 0.2"
end select
end do
end program anime_birdseyeview_rotation_generator
上に述べた市松模様状に+4と-4の加熱と冷却の分布を与えたときの熱平衡温度分布を求めるプログラムを thermal_diffusion.f90を変更して作り、それをthermal_diffusion_checker.f90とせよ。
as of 2010-08-02 (月) 16:11:15 (348)