- 追加された行はこの色です。
- 削除された行はこの色です。
-------
担当： 政田 洋平
-------
【目次】
#contents
-------
* 可視化の必要性 [#pfc0c233]
人間の視覚を使って情報を引き出す操作を一般に可視化という。
科学技術計算においては、可視化が様々な局面で重要な役割を果たす。
特に、計算結果として出力される数値データに潜む情報を引き出すための可視化―これをデータ可視化という―は、
不可欠と言ってよい。
データ可視化の身近な例は天気図であろう。
NKHラジオの「気象通報」という番組を知っているだろうか？
各地の風向、風力、気圧などの数字をアナウンサーが次々に読み上げていくラジオ番組である。
これをただ耳で聞いているだけでは何もわからない。
しかしペンと地図を手にし、番組を聞きながら読み上げられる各地点に気圧などを書き込み（専用の地図を売っている）、
最後に同じ気圧になりそうな点を曲線（等高線）でつなげていくと、
低気圧や高気圧が図―天気図―としてあらわれる。
知識のある人が天気図を見れば、 明日の天気がだいたいわかる。
この天気図に入っている情報は、
もともとラジオのアナウンサーが読み上げる数値データに全て含まれているものである。
だが、そのような数値データの羅列をいくら注意深く聞いても、
あるいはそれを数表にして紙の上に書き、それをじっと睨んでみても、
明日の天気を知るのはほとんど不可能である。
紙の上に天気図という図にしてみてはじめてデータに隠されていた情報を引き出すことができる。
これが可視化である。
* 可視化の例 [#h364e87b]
データの可視化には可視化の目的や、対象データの種類の応じて様々な方法が開発されている。
その中で、ここでは、スカラー場の可視化について簡単に説明する。
** 1次元グラフ [#za71e69f]
xのある関数f(x)を直感的に理解するためには、グラフを書くのが一番である。
例えばxのx乗の関数、つまりf(x)=x^xはどんな形だろうか？
下のグラフを見ればわかる。
#ref(graph_1d.png)
** 2次元等高線 [#t928cb5d]
xとyの関数、つまり2次元の関数f(x,y)の形を理解するには、等高線を描くのがよい。
地表面での大気の圧力pの分布p(x,y)の等高線は天気図でおなじみである。
#ref(graph_2d.png)
上の図では等高線の高さを色で表示してさらに分かりやすくしている。
ここで等高線の描画アルゴリズムを紹介しよう。
計算格子点上に定義されたデータから一本の等高線を描くには、下の図のように短い線分をつなげていけばよい。
#ref(contour_line01.png)
一つの線分は4つの計算格子点で定義された長方形領域（セル）の中で直線を描く。
例えばf(x,y)=1.0の値の等高線を描く場合を考えよう。
あるセルの4つの頂点におけるfの値が全て1.0よりも大きいか、あるいは1.0未満であれば、f=1.0の等高線はこのセルを通らないのは明らかである。
セルを周囲の4つの辺、それぞれの両端の頂点でのfの値が1.0を「挟めば」その辺を等高線が通る。
辺上のどの位置を等高線が横切るかは、線形補間をすればよい。
下の図はちょうど中点を通る例である。
#ref(contour_line02.png)
このアルゴリズムはmarching squaresと呼ばれる。
等高線による可視化とは2次平面上に分布するスカラー場f(x,y)を
等高線という曲線を見ることで理解する方法である。
** 3次元等値面 [#va1d6f14]
等値面の「3次元版」も考えられる。
この場合は3次元空間中の曲面の形状を見ることで3次元空間に分布する関数f(x,y,z)の分布を理解する。
たとえば、f(x,y,z)がある値（例えば1.0）をとる点の集合は、曲面の方程式 f(x,y,z)=1.0 で定義される。
下の図の白い矢印で示している物体は等値面の例である。
これは核融合実験装置（LHD）内部のプラズマの圧力分布の等値面である。
#ref(isosurface_sample.png)
このような等値面を描くためのアルゴリズムとしてMarching Cubesがある。
#ref(marching_cubes.png)
このMarching Cubesアルゴリズムは、
1985年に米国で特許となったが、20年を経過したので、現在では自由に利用できる。
* gnuplot入門 [#zd84152f]
このような可視化アルゴリズムを実装したソフトウェアは既に各種開発されており、可視化ソフトウェアとかグラフソフトなどと呼ばれる。
この演習では、そのなかで無償で利用出来る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で計算した正方形領域の中心点での温度が、
ヤコビ法によって次第に収束して行く様子を可視化してみよう。
* &color(#0000ff){【演習】1次元グラフ１}; [#pb2617b4]
(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と書いてもよい。）
** 結果例 [#nd310634]
#ref(central_temp_converge.png)
* &color(#0000ff){【演習】1次元グラフ２}; [#zb0c8a9d]
(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
と打て。
** 結果例 [#e2532ba8]
#ref(central_temp_converge_rev.png)
* &color(#0000ff){【演習】1次元グラフ３}; [#fdc3f653]
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という名前のファイルが出来ているはずである。
* &color(#0000ff){【演習】1次元グラフ４}; [#q7d70abf]
../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
と打て。
** 結果例 [#qe50d43a]
#ref(temp_final_x_profile.png)
(4) ファイルheat1_print_final_x_prof.gpを自由に変更してグラフがどう変わるか試せ。
なお、gnuplotのヘルプはコマンドプロンプトでhelpと打てば得られる。
* アニメーション用データの生成 [#u4f89dc0]
次にアニメーションによって収束の様子を確認しよう。
そのためのデータ（連番つきファイル群）を書き出すためのプログラム
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によるアニメーション [#o7126c27]
gnuplotのスクリプトを使えばアニメーションも簡単にできる。
ただし、コマ数が多くなると手でスクリプトを書くのは難しくなるので、
「gnuplotスクリプト自動生成プログラム」を作ったほうがよい。
以下はその例である。
* &color(#0000ff){【演習】1次元グラフ アニメーション}; [#b2b4e3cd]
以下のコードを
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
* 先週の熱伝導問題を解くプログラムの復習 [#f6966930]
- コードの構造（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次元可視化 [#k39ed7ce]
以降では、2次元データの可視化を行う。これまでと同様gnuplotを利用する。
gnuplotで2次元データをプロットする場合、そのデータのフォーマットは
x00 y00 関数値
x01 y00 関数値
x02 y00 関数値
・
・
・
x09 y00 関数値
（空行）
x00 y01 関数値
x01 y01 関数値
x02 y01 関数値
・
・
・
x09 y01 関数値
（空行）
・
・
・
x09 y09 関数値
である。
#br
#br
#br
* 2次元データ出力ルーチン [#a2717789]
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
#br
#br
#br
** 出力データの確認 [#m11f0188]
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
・
・
・
#br
#br
#br
* &color(#0000ff){【演習】2次元等高線の表示}; [#n213fe1f]
- 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
- 結果の例
#ref(thermal_diffusion_2d_contour_lines.png)
#br
#br
#br
* &color(#0000ff){【演習】アニメーション}; [#n213fe1f]
- 時間発展（アニメーション）で等高線を見てみよう。
- 例によってスクリプト生成プログラムを使う。
- ファイル名：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
#br
#br
#br
* &color(#0000ff){【演習】色分布による可視化（静止画）}; [#n213fe1f]
- 等高線を描く代わりに正方形領域内部各点の温度を色で表現することも可能である。
- 実際に描いてみよう。
- 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
- 結果の例
#ref(thermal_diffusion_2d_contour_colors.png)
#br
#br
#br
* &color(#0000ff){【演習】色分布による可視化（アニメーション）}; [#n213fe1f]
- 時間発展のアニメーションを作ってみよう。
- スクリプト生成プログラム 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
#br
#br
#br
* &color(#0000ff){【演習】gnuplotによる鳥瞰図（静止画）}; [#n213fe1f]
- 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
- 結果の例
#ref(thermal_diffusion_2d_birdseyeview.png)
#br
#br
#br
* &color(#0000ff){【演習】gnuplotによる鳥瞰図（時間発展アニメーション）}; [#n213fe1f]
- アニメーションを作ってみよう。
- ファイル名：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
#br
#br
#br
* &color(#0000ff){【演習】gnuplotによる鳥瞰図（回転のアニメーション）}; [#n213fe1f]
- 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
#br
#br
#br
* 応用問題： 熱源分布が一様でない場合の平衡温度分布 [#w71b98ed]
- これまで正方形領域（0<=x<=1, 0<=y<=1）が全体に一様に+4の強さで加熱されている場合を想定して計算していた。
#ref(heating_uniform.png)
- 次に、下の図のように正方形を4つの部分にわけて市松模様状に加熱と冷却の分布を与えることを考える。
冷却は負の加熱と考え、+4の加熱と-4の加熱の分布を図のように与える。
#ref(heating_checker.png)
- 答えの分布をsample_birdseyeview.gpで可視化すると次のような図が得られる。
#ref(thermal_diffusion_checker_birdseyeview_sample_stamp.png)
#br
#br
#br
* &color(#ff0000){レポート課題}; [#od9affb6]
【提出は次週までにメールで（宛先：ymasada_at_harbor.kobe-u.ac.jp）】
~
熱伝導問題を解くコードを応用し、上に述べた市松模様状に+4と-4の加熱と冷却の分布を与えたときの熱平衡温度分布を求めるプログラムを作成せよ。
+ 市松模様状の加熱と冷却を与えるために、熱伝導問題を解くコードにどのような変更を加えたか説明せよ。
+ 市松模様状に加熱と冷却を与えた際の平衡解（最終状態解）温度の最大値（正の値）と最小値（負の値）を記せ。
+ その平衡解の温度分布を可視化した図をつけること。gnuplotでなくても可。なおファイル形式はpdf, png, jpgのいずれか。ファイルサイズは100kB以下にすること。
* 授業アンケート [#veca5581]
（どれか一つ、一度だけ押してください。）
難易度
#vote(簡単すぎた[7], 難しすぎた[12], ちょうどよかった[10])
分量
#vote(少ない[1], ちょうどよい[4], 少し多い[9], 多すぎる[1])
可視化ソフトについて
#vote(gnuplotでよかった[11], 別のグラフソフトの方がよい[3], Processingを使うべき[1], なんでもよい[6])
#vote(gnuplotでよかった[11], 別のグラフソフトの方がよい[3], Processingを使うべき[1], なんでもよい[7])
* 質問、コメントなど自由にどうぞ [#q4a8df41]
「お名前」欄は空欄で可。
#pcomment
------------------------------
as of &_now; (&counter;)
