6.1 1次元可視化 のバックアップ(No.2)
- バックアップ一覧
- 差分 を表示
- 現在との差分 を表示
- ソース を表示
- 6.1 1次元可視化 へ行く。
- 1 (2012-06-28 (木) 14:15:26)
- 2 (2012-06-28 (木) 14:19:56)
現在（2013-07-02 (火) 18:18:57）作成中です。 既に書いてある内容も大幅に変わる可能性が高いので注意。
神戸大学 大学院システム情報学研究科 計算科学専攻 陰山 聡
【目次】
可視化の必要性 †
上の例では、MPIプロセスがもつ配列uの値をテキストデータとして書き出すことで、 プロセス間のMPI_SENDRECVによる通信の様子が分かりやすく表示された。 このように人間の視覚を使って情報を引き出す操作を一般に可視化という。
科学技術計算においては、可視化が様々な局面で重要な役割を果たす。 特に、計算結果として出力される数値データに潜む情報を引き出すための可視化―これをデータ可視化という―は、 不可欠と言ってよい。
データ可視化の身近な例は天気図であろう。
NKHラジオの「気象通報」という番組を知っているだろうか？ 各地の風向、風力、気圧などの数字をアナウンサーが次々に読み上げていくラジオ番組である。 これをただ耳で聞いているだけでは何もわからない。 しかしペンと地図を手にし、番組を聞きながら読み上げられる各地点に気圧などを書き込み（専用の地図を売っている）、 最後に同じ気圧になりそうな点を曲線（等高線）でつなげていくと、 低気圧や高気圧が図―天気図―としてあらわれる。 知識のある人が天気図を見れば、 明日の天気がだいたいわかる。
この天気図に入っている情報は、 もともとラジオのアナウンサーが読み上げる数値データに全て含まれているものである。 だが、そのような数値データの羅列をいくら注意深く聞いても、 あるいはそれを数表にして紙の上に書き、それをじっと睨んでみても、 明日の天気を知るのはほとんど不可能である。 紙の上に天気図という図にしてみてはじめてデータに隠されていた情報を引き出すことができる。 これが可視化である。
可視化の例 †
データの可視化には可視化の目的や、対象データの種類の応じて様々な方法が開発されている。 その中で、ここでは、スカラー場の可視化について簡単に説明する。
1次元グラフ †
xのある関数f(x)を直感的に理解するためには、グラフを書くのが一番である。 例えばxのx乗の関数、つまりf(x)=x^xはどんな形だろうか？ 下のグラフを見ればわかる。
2次元等高線 †
xとyの関数、つまり2次元の関数f(x,y)の形を理解するには、等高線を描くのがよい。 地表面での大気の圧力pの分布p(x,y)の等高線は天気図でおなじみである。
上の図では等高線の高さを色で表示してさらに分かりやすくしている。
ここで等高線の描画アルゴリズムを紹介しよう。 計算格子点上に定義されたデータから一本の等高線を描くには、下の図のように短い線分をつなげていけばよい。
一つの線分は4つの計算格子点で定義された長方形領域（セル）の中で直線を描く。 例えばf(x,y)=1.0の値の等高線を描く場合を考えよう。 あるセルの4つの頂点におけるfの値が全て1.0よりも大きいか、あるいは1.0未満であれば、f=1.0の等高線はこのセルを通らないのは明らかである。
セルを周囲の4つの辺、それぞれの両端の頂点でのfの値が1.0を「挟めば」その辺を等高線が通る。 辺上のどの位置を等高線が横切るかは、線形補間をすればよい。 下の図はちょうど中点を通る例である。
このアルゴリズムは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）内部のプラズマの圧力分布の等値面である。
このような等値面を描くためのアルゴリズムとしてMarching Cubesがある。
この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と打つ。
gnuplotではこのように解析的に定義された関数をターミナルに打ち込んでそのグラフを描く機能がある。 それに加えて、ファイルに書き込まれた離散データを読み込み、それをグラフにする機能もある。 この機能を使ってheat1.f90で計算した正方形領域の中心点での温度が、 ヤコビ法によって次第に収束して行く様子を可視化してみよう。
【演習】1次元グラフ１ †
(1) heat1.f90の出力（標準出力）をファイルtest.dataという名前に保存せよ。
- コンパイルはpgf95 heat1.f90
- ./a.out > test.data でリダイレクトされる。
(2) test.dataの中身を確認せよ
- エディタで開くよりもheadコマンドやtailコマンドで見る方が早い。
(3) gnuplotを立ち上げ、コマンドプロンプトに plot 'test.data' w lp と入れよ。 lpはlinespointsの略で、線（line）と点（point）を表示することを意味する。 （w linespointsと書いてもよい。）
結果例 †
【演習】1次元グラフ２ †
(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
と打て。
結果例 †
【演習】1次元グラフ３ †
heat1.f90で計算された最終的な平衡温度分布をgnuplotで見てみよう。 正方形を真ん中で横に切るy=0.5の線上での温度のx分布をグラフにする。
(1) これからの演習でデータファイルを多数生成するので、データファイル出力先のディレクトリを用意しよう。 例えばいまソースコードの置いていあるディレクトリと同じ高さに
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))/4.0_DP+heat*h*h
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次元グラフ４ †
../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
と打て。
結果例 †
(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))/4.0_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
レポート課題 †
・・・未定・・・
授業アンケート †
今回の演習内容はどうでしたか？（どれか一つ、一度だけ押してください。）
質問、コメントなど自由にどうぞ †
「お名前」欄は空欄で可。
コメントはありません。 コメント/6.1 1次元可視化?
as of 2023-09-23 (土) 00:11:53 (1420)