担当: 政田 洋平


【目次】


可視化の必要性

人間の視覚を使って情報を引き出す操作を一般に可視化という。

科学技術計算においては、可視化が様々な局面で重要な役割を果たす。 特に、計算結果として出力される数値データに潜む情報を引き出すための可視化―これをデータ可視化という―は、 不可欠と言ってよい。

データ可視化の身近な例は天気図であろう。

NKHラジオの「気象通報」という番組を知っているだろうか? 各地の風向、風力、気圧などの数字をアナウンサーが次々に読み上げていくラジオ番組である。 これをただ耳で聞いているだけでは何もわからない。 しかしペンと地図を手にし、番組を聞きながら読み上げられる各地点に気圧などを書き込み(専用の地図を売っている)、 最後に同じ気圧になりそうな点を曲線(等高線)でつなげていくと、 低気圧や高気圧が図―天気図―としてあらわれる。 知識のある人が天気図を見れば、 明日の天気がだいたいわかる。

この天気図に入っている情報は、 もともとラジオのアナウンサーが読み上げる数値データに全て含まれているものである。 だが、そのような数値データの羅列をいくら注意深く聞いても、 あるいはそれを数表にして紙の上に書き、それをじっと睨んでみても、 明日の天気を知るのはほとんど不可能である。 紙の上に天気図という図にしてみてはじめてデータに隠されていた情報を引き出すことができる。 これが可視化である。

可視化の例

データの可視化には可視化の目的や、対象データの種類の応じて様々な方法が開発されている。 その中で、ここでは、スカラー場の可視化について簡単に説明する。

1次元グラフ

xのある関数f(x)を直感的に理解するためには、グラフを書くのが一番である。 例えばxのx乗の関数、つまりf(x)=x^xはどんな形だろうか? 下のグラフを見ればわかる。

graph_1d.png

2次元等高線

xとyの関数、つまり2次元の関数f(x,y)の形を理解するには、等高線を描くのがよい。 地表面での大気の圧力pの分布p(x,y)の等高線は天気図でおなじみである。

graph_2d.png

上の図では等高線の高さを色で表示してさらに分かりやすくしている。

ここで等高線の描画アルゴリズムを紹介しよう。 計算格子点上に定義されたデータから一本の等高線を描くには、下の図のように短い線分をつなげていけばよい。

contour_line01.png

一つの線分は4つの計算格子点で定義された長方形領域(セル)の中で直線を描く。 例えばf(x,y)=1.0の値の等高線を描く場合を考えよう。 あるセルの4つの頂点におけるfの値が全て1.0よりも大きいか、あるいは1.0未満であれば、f=1.0の等高線はこのセルを通らないのは明らかである。

セルを周囲の4つの辺、それぞれの両端の頂点でのfの値が1.0を「挟めば」その辺を等高線が通る。 辺上のどの位置を等高線が横切るかは、線形補間をすればよい。 下の図はちょうど中点を通る例である。

contour_line02.png

このアルゴリズムはmarching squaresと呼ばれる。 等高線による可視化とは2次平面上に分布するスカラー場f(x,y)を 等高線という曲線を見ることで理解する方法である。

3次元等値面

等値面の「3次元版」も考えられる。 この場合は3次元空間中の曲面の形状を見ることで3次元空間に分布する関数f(x,y,z)の分布を理解する。 たとえば、f(x,y,z)がある値(例えば1.0)をとる点の集合は、曲面の方程式 f(x,y,z)=1.0 で定義される。

下の図の白い矢印で示している物体は等値面の例である。 これは核融合実験装置(LHD)内部のプラズマの圧力分布の等値面である。

isosurface_sample.png

このような等値面を描くためのアルゴリズムとしてMarching Cubesがある。

marching_cubes.png

このMarching Cubesアルゴリズムは、 1985年に米国で特許となったが、20年を経過したので、現在では自由に利用できる。

gnuplot入門

このような可視化アルゴリズムを実装したソフトウェアは既に各種開発されており、可視化ソフトウェアとかグラフソフトなどと呼ばれる。 この演習では、そのなかで無償で利用出来るgnuplotを利用する。 (なお、gnuplotはGNUのソフトウェアではない。)

/usr/bin/gnuplot

または

/mnt/nfs/usr/bin/gnuplot

というコマンド打てばgnuplotが立ち上がる。

gnuplotコマンドプロンプトに

plot sin(x)

と打ち込んでみよう。sin関数がプロットされるはずである。 上に例として示したxのx乗のグラフは plot x**xとすれば見える。

gnuplotの終了はコマンドプロンプトでquit(またはexit)と打つ。

gnuplotではこのように解析的に定義された関数をターミナルに打ち込んでそのグラフを描く機能がある。 それに加えて、ファイルに書き込まれた離散データを読み込み、それをグラフにする機能もある。 この機能を使ってtemperature.f90で計算した正方形領域の中心点での温度が、 ヤコビ法によって次第に収束して行く様子を可視化してみよう。

【演習】1次元グラフ1

(1) heat1c.f90の出力(標準出力)をファイルtest.dataという名前に保存せよ。

  • コンパイルはpgf95 heat1c.f90
  • ./a.out > test.data でリダイレクトされる。

(2) test.dataの中身を確認せよ

  • エディタで開くよりもheadコマンドやtailコマンドで見る方が早い。

(3) gnuplotを立ち上げ、コマンドプロンプトに plot 'test.data' w lp と入れよ。 lpはlinespointsの略で、線(line)と点(point)を表示することを意味する。 (w linespointsと書いてもよい。)

結果例

central_temp_converge.png

【演習】1次元グラフ2

(1) 上の図の中で文字とグラフが重なって見にくいので、文字を消す。 gnuplotのコマンドプロンプトで

unset key

と入れ、それに続いて

replot

と打て。

(2) 縦軸の表示範囲を調整をする。 gnuplotのコマンドプロンプトで

set yrange [0:0.35]
replot

と打て。

(3) 図全体のタイトルと、x軸、y軸の説明を入れる。

set title "Temperatute at the Center"
set xlabel "iterations"
set ylabel "temperature"
replot

と打て。

結果例

central_temp_converge_rev.png

【演習】1次元グラフ3

heat1c.f90で計算された最終的な平衡温度分布をgnuplotで見てみよう。 正方形を真ん中で横に切るy=0.5の線上での温度のx分布をグラフにする。

(1) これからの演習でデータファイルを多数生成するので、データファイル出力先のディレクトリを用意しよう。 例えばいまソースコード(heat1c.f90)の置いていあるディレクトリと同じ高さで

mkdir ../data

としてdataディレクトリを作成する。別の場所、名前でも構わないが、以後のサンプルコードではそれに応じて適宜dataディレクトリの名前を変更すること。

(2) 以下のプログラムをheat1_print_final_x_prof.f90 という名前のソースコードにコピーし、コンパイル・実行せよ。

program heat1_print_final_x_prof
  implicit none
  integer, parameter :: m=31, nmax=5000
  integer :: i,j,n
  integer, parameter :: SP = kind(1.0)
  integer, parameter :: DP = selected_real_kind(2*precision(1.0_SP))
  real(DP), dimension(:,:), allocatable :: u, un
  real(DP) :: h, heat=1.0_DP, heat_hsq
  allocate(u(0:m+1,0:m+1))
  allocate(un(m,m))
  h = 1.0_DP/(m+1)
  heat_hsq = heat*h*h
  u(:,:) = 0.0_DP
  open(10,file='../data/temp.final_profile_x')
  do n = 1 , nmax
    do j = 1 , m
      do i = 1 , m
        un(i,j)=(u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1))*0.25_DP+heat_hsq
      end do
    end do
    u(1:m,1:m)=un(1:m,1:m)
  end do
  do i = 0 , m+1
     write(10,*) i, u(i,m/2+1)
  end do
  deallocate(u,un)
  close(10)
end program heat1_print_final_x_prof

ディレクトリ../dataにtemp.final_profile_xという名前のファイルが出来ているはずである。

【演習】1次元グラフ4

../data/temp.final_profile_xのデータをgnuplotでグラフにする。

(1) このデータの中身をエディタまたはheadコマンド、tailコマンドで確認せよ。

(2) gnuplotではコマンドプロンプトに入力する内容をスクリプトファイルから読み込ませることが出来る。 以下の内容をheat1_print_final_x_prof.gpという名前のファイルに保存せよ。

#
# heat1_print_final_x_prof.gp
#
# final temperature profile at y=0.5 as a function of x
#
set xrange [0:32]
set yrange [0:0.5]
set xlabel 'i'
set ylabel 'temp at y=0.5'
plot '../data/temp.final_profile_x' w lp 
pause -1

(3) gnuplotがまだ立ち上がっていたらquitコマンドで終了し、改めてshellで

gnuplot heat1_print_final_x_prof.gp

と打て。

結果例

temp_final_x_profile.png

(4) ファイルheat1_print_final_x_prof.gpを自由に変更してグラフがどう変わるか試せ。 なお、gnuplotのヘルプはコマンドプロンプトでhelpと打てば得られる。

アニメーション用データの生成

次にアニメーションによって収束の様子を確認しよう。 そのためのデータ(連番つきファイル群)を書き出すためのプログラム heat_animation_x_prof.f90 を以下に示す。 このプログラムをpgf90コマンドでコンパイルし、 実行すればdataディレクトリに連番ファイルが出力される。

module constants
  implicit none
  integer, parameter :: m=61, nmax=5000
  integer, parameter :: SP = kind(1.0)
  integer, parameter :: DP = selected_real_kind(2*precision(1.0_SP))
end module constants
module print
  use constants
  implicit none
contains
  subroutine print__profile(f)
    real(DP), dimension(0:m+1), intent(in) :: f
    integer :: i
    integer :: counter = 0  ! automatic save attribute
    character(len=*), parameter :: base = "../data/temp.j=middle."
    character(len=4) :: serial_num
    write(serial_num,'(i4.4)') counter
    open(10,file=base//serial_num)
    do i = 0 , m+1
       write(10,*) i, f(i)
    end do
    close(10)
    counter = counter + 1
  end subroutine print__profile
end module print
program heat1_animation_x_prof
  use constants
  use print
  implicit none
  integer :: i,j,n
  real(DP), dimension(:,:), allocatable :: u, un
  real(DP) :: h, heat=1.0_DP
  real(DP) :: random
  allocate(u(0:m+1,0:m+1))
  allocate(un(m,m))
  h = 1.0_DP/(m+1)
  u(:,:) = 0.0_DP
  call print__profile(u(:,m/2+1))
  do n=1, nmax
    do j=1, m
      do i=1, m
        un(i,j)=(u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1))*0.25_DP+heat*h*h
      end do
    end do
    u(1:m,1:m)=un(1:m,1:m)
    if (mod(n,100)==0) call print__profile(u(:,m/2+1))
  end do
  deallocate(u,un)
end program heat1_animation_x_prof

gnuplotによるアニメーション

gnuplotのスクリプトを使えばアニメーションも簡単にできる。 ただし、コマ数が多くなると手でスクリプトを書くのは難しくなるので、 「gnuplotスクリプト自動生成プログラム」を作ったほうがよい。 以下はその例である。

【演習】1次元グラフ アニメーション

以下のコードを heat1_animation_x_prof_gp_generator.f90 というファイルに保存し、ソースコードの冒頭にあるUsageに従ってgnuplotによる可視化アニメーションを作れ。

!
! heat1_animation_x_prof_gp_generator.f90
!
!    by  Akira Kageyama
!    at  Kobe University
!    on  2010.07.07
!   for  data visualization of heat1_animation_x_prof.f90
!    in  the Lecture "Computational Science" (2010)
!
! Usage: 
!    (1) check the mesh size m and the counter_end.
!    (2) pgf95 this_code
!    (3) ./a.out > anyname
!    (4) gnuplot anyname
!

program heat1_animation_x_prof_gp_generator
  implicit none
  integer, parameter :: m = 61
  integer :: counter, counter_end = 50
  character(len=*), parameter :: base='../data/temp.j=middle.'
  character(len=4)            :: serial_num
  
  print *, "#"
  print *, "# gnuplot script generated by heat1_animation_x_prof_gp_generator.f90"
  print *, "#"
  print *, " "
  print *, "set xlabel 'j'              # x-axis"
  print *, "set ylabel 'temperature'    # y-axis"   
  print *, "set xrange[0:", m+1, "]     # j-grid min & max"
  print *, "set yrange[0.0:0.5]         # temperature min & max"

  do counter = 0 , counter_end
     write(serial_num,'(i4.4)') counter
     print *, "plot '"//base//serial_num//"' w lp"
     if ( counter==0) then 
        print *, "pause 5"
     else
        print *, "pause 1"
     end if
  end do
  print *, "pause -1"
end program heat1_animation_x_prof_gp_generator

先週の熱伝導問題を解くプログラムの復習

  • コードの構造(120705/heat7 ディレクトリを見よ)
  • モジュール構成は以下のとおり:
モジュール名内容
ut (constants含む)便利ツール
stopwatch計算時間測定
parallel並列化
temperature温度場
  • main programの構造

データ出力、時間計測、assertのルーチンを除くとコードの構造は以下のとおり:

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

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

  • 先週説明した熱伝導問題を解くコードでは、例えばmainプログラム
 call ut__assert()
 call parallel__initialize
 call temperature__initialize
 call stopwatch__stt()

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

2次元可視化

以降では、2次元データの可視化を行う。これまでと同様gnuplotを利用する。

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

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

である。

 
 
 

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

temperature.f90にある2次元データ出力ルーチンoutput_2d_profile は、正方形上(x,y平面上)に分布する温度をすべて書き出すものである。 中身を見てみよう。

 subroutine output_2d_profile
   real(DP), dimension(0:NGRID+1,    &
                       0:NGRID+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 ( p%myrank==0 ) then
      open(10,file=base//serial_num)
      do j = 0 , NGRID+1
         do i = 0 , NGRID+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 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 (heat7/visに置いてある)
  • 実行方法: 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_sample_stamp.png
 
 
 

レポート課題

【提出は次週までにメールで(宛先:ymasada_at_harbor.kobe-u.ac.jp)】
熱伝導問題を解くコードを応用し、上に述べた市松模様状に+4と-4の加熱と冷却の分布を与えたときの熱平衡温度分布を求めるプログラムを作成せよ。

  1. 市松模様状の加熱と冷却を与えるために、熱伝導問題を解くコードにどのような変更を加えたか説明せよ。
  2. 市松模様状に加熱と冷却を与えた際の平衡解(最終状態解)温度の最大値(正の値)と最小値(負の値)を記せ。
  3. その平衡解の温度分布を可視化した図をつけること。gnuplotでなくても可。なおファイル形式はpdf, png, jpgのいずれか。ファイルサイズは100kB以下にすること。

授業アンケート

(どれか一つ、一度だけ押してください。)

難易度

選択肢 投票
簡単すぎた 9  
難しすぎた 13  
ちょうどよかった 12  

分量

選択肢 投票
少ない 1  
ちょうどよい 4  
少し多い 9  
多すぎる 2  

可視化ソフトについて

選択肢 投票
gnuplotでよかった 13  
別のグラフソフトの方がよい 3  
Processingを使うべき 2  
なんでもよい 7  

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

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

コメントはありません。 コメント/7.データ可視化?

お名前:

as of 2018-07-20 (金) 23:03:02 (12939)