最終更新:2019-10-15 (火) 01:08:25


神戸大学 大学院システム情報学研究科 計算科学専攻 陰山 聡

 
  • 先週は私の海外出張のため、政田先生に代わっていただいた。
  • 今週と来週のあと2回。
  • 今回使うソースコードは、scalar の /tmp/110714に置いた。
  • wiki中、ソースコード部分が等角フォントで表示されていない場合、 別のブラウザを試してみること。 (ソースコード中のコメントは一部、等角フォントを仮定している。)

【目次】


 
 
 

問題設定の復習

  • 2次元正方形領域 [0,1]×[0,1] での熱伝導を考える
  • 境界をすべて0℃に固定
  • 領域全体に一定の熱を加える
  • ⇒ このとき,十分な時間が経った後での温度分布はどうなるか?
 
 
 

前回の宿題

前回の宿題は、「thermal_diffusion.f90を理解しておくこと」であった。 ここではポイントだけ説明する。

  • コード全体の構造
  • モジュール構成は以下のとおり:
モジュール名内容
constants定数
parallel並列化
temperature温度場
  • main programの構造

データ出力ルーチンを除くとコードの構造は以下のとおり:

 call parallel__initialize
 call temperature__initialize
 do loop = 1 , LOOP_MAX
    call temperature__update
 end do
 call temperature__finalize
 call parallel__finalize

ルーチン名を見れば意味は明白であろう。

  • thermal_diffusion.f90コードでは、例えばmainプログラム
 call parallel__initialize
 call temperature__initialize
 call temperature__output_1d_profile
 call temperature__output_2d_profile

に見られるように、関数名(サブルーチン名)は、「モジュール名」+「アンダースコア二つ」+「その関数の中身を示す言葉」で統一されている。 その関数がどのモジュールで定義されているかがすぐに分かるようにするために工夫である。

  • 加熱項の表式

今考えている正方形領域は 一様かつ一定に加熱されている。 その加熱率 heat_source を(ある規格化の下で)4とすると、 定常状態における温度T(x,y)は以下のポアッソン方程式を満たす。

poisson_eq_01.png
 
 
 
  • ちょっと確認: ラプラシアン記号を使った上の式の意味は分かりますか?
    選択肢 投票
    もちろん分かる 2  
    ちょっと不安がある 8  
    全く知らない 11  

上の式を丁寧に書くと、

poisson_eq_02.png

である。

 
 
 
  • もう一つ確認: 「何か」のラプラシアンをとると、その「 何か」のガタガタをならす効果があります。 この意味が分かりますか?
    選択肢 投票
    もちろん分かる 1  
    ちょっと不安がある 3  
    さっぱり分からない 13  

連続場としての温度場を 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を繰り返し計算する。

  • heat_sourceが一定かつ一様で4という値であれば、右辺第二項は
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という定数を足すだけで済むようにしている。

  • 上のコードをよく見ると、コメントアウトされた部分では4.0で割っているところを、 コメントアウトされていない行では0.25を掛けている。 これは割り算よりも掛け算の方が計算が速いからである。
 
 
 

【演習】 テストデータの作成

  • thermal_diffusion.f90をmpif90コマンドでコンパイルし、
  • ジョブスクリプトthermal_diffusion.jsで並列計算せよ。
 
 
 

2次元可視化

今回は2次元データの可視化を行う。 前回と同様gnuplotを利用する。

gnuplotで2次元データをプロットする場合、そのデータのフォーマットは

x00 y00 関数値
x01 y00 関数値
x02 y00 関数値
 ・
 ・
 ・
x09 y00 関数値
(空行)
x00 y01 関数値
x01 y01 関数値
x02 y01 関数値
 ・
 ・
 ・
x09 y01 関数値
(空行)
 ・
 ・
 ・
x09 y09 関数値

である。

 
 
 

2次元データ出力ルーチン

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     
          ・
          ・
          ・
 
 
 

【演習】2次元等高線の表示

  • data/temp.2d.0100のファイルに記された温度の分布をgnuplotの等高線で可視化してみよう。
  • ファイル名:sample_contour_lines.gp
  • 実行方法: gnuplot sample_contour_lines.gp
  • ファイル名やパラメータ等を自由に変更してその効果を試せ。
#
# 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

  • 結果の例
    thermal_diffusion_2d_contour_lines.png
 
 
 

【演習】アニメーション

  • 時間発展(アニメーション)で等高線を見てみよう。
  • 例によってスクリプト生成プログラムを使う。
  • ファイル名:anime_contour_lines_generator.f90
  • 利用方法は以下のコメント(Usage)に書かれている。
!
! 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
 
 
 

【演習】色分布による可視化(静止画)

  • 等高線を描く代わりに正方形領域内部各点の温度を色で表現することも可能である。
  • 実際に描いてみよう。
  • gnuplotのサンプルスクリプトは以下のとおり。ファイル名: sample_contour_colors.gp
#
# 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
  • 結果の例
    thermal_diffusion_2d_contour_colors.png
 
 
 

【演習】色分布による可視化(アニメーション)

  • 時間発展のアニメーションを作ってみよう。
  • スクリプト生成プログラム anime_contour_colors_generator.f90 (利用方法はUsage以下の通り)
!
! 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
 
 
 

【演習】gnuplotによる鳥瞰図(静止画)

  • 2次元温度分布T(x,y)を高さ(z)で表すことも可能である。 この時の描画は3次元的に行う必要がある。 空をとぶ鳥から見下ろしたような図は一般に鳥瞰図(bird's eye view)とも呼ばれる。
  • gnuplotで描いてみよう。
  • ファイル名:sample_birdseyeview.gp
#
# 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

  • 結果の例
    thermal_diffusion_2d_birdseyeview.png
 
 
 

【演習】gnuplotによる鳥瞰図(時間発展アニメーション)

  • アニメーションを作ってみよう。
  • ファイル名:anime_birdseyeview_generator.f90
!
! 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
 
 
 

【演習】gnuplotによる鳥瞰図(回転のアニメーション)

  • gnuplotではviewというパラメータで視線の方向(「鳥」がどの方向から見下ろしているか)を指定することが出来る。 このviewパラメータを変更したアニメーションを作ってみよう。
  • ファイル名:anime_birdseyeview_rotation_generator.f90
!
! 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
 
 
 

応用問題: 熱源分布が一様でない場合の平衡温度分布

  • これまで正方形領域(0<=x<=1, 0<=y<=1)が全体に一様に+4の強さで加熱されている場合を想定して計算していた。
    heating_uniform.png
  • 次に、下の図のように正方形を4つの部分にわけて市松模様状に加熱と冷却の分布を与えることを考える。 冷却は負の加熱と考え、+4の加熱と-4の加熱の分布を図のように与える。
    heating_checker.png
  • 答えの分布をsample_birdseyeview.gpで可視化すると次のような図が得られる。
thermal_diffusion_checker_birdseyeview.png
 
 
 

レポート課題 1

上に述べた市松模様状に+4と-4の加熱と冷却の分布を与えたときの熱平衡温度分布を求めるプログラムを thermal_diffusion.f90を変更して作り、それをthermal_diffusion_checker.f90とせよ。

  1. thermal_diffusion_checker.f90を作るにあたり、thermal_diffusion.f90にどのような変更を加えたか説明せよ。
  2. thermal_diffusion_checker.f90で求めた平衡解(最終状態解)温度の最大値と最小値を記せ。
  3. その最大値と最小値が格子分割点数 MESH_SIZEにどう依存性するか調べ、考察せよ。
  4. 【来週の講義で提出】
 
 
 

授業アンケート

今回の演習内容はどうでしたか?(どれか一つ、一度だけ押してください。)

選択肢 投票
簡単すぎた 0  
難しすぎた 3  
ちょうどよかった 10  
 
 
 

質問、コメントなど自由にどうぞ

「お名前」欄は空欄で可。

最新の10件を表示しています。 コメントページを参照

  • 感想、質問など気楽に書き込んでください。もちろん匿名で可。名前欄を空白のままでOK。 -- 陰山? 2010-07-15 (木) 00:13:47
  • これはテストです。 -- テスト? 2011-03-11 (金) 14:35:26
お名前:

as of 2019-10-15 (火) 01:08:25 (5604)