[#pfc0c233]人間の視覚を使って情報を引き出す操作を一般に可視化という。科学技術計算においては、可視化が様々な局面で重要な役割を果たす。データ可視化の身近な例は天気図であろう。NKHラジオの「気象通報」という番組を知っているだろうか?各地の風向、風力、気圧などの数字をアナウンサーが次々に読み上げていくラジオ番組である。これをただ耳で聞いているだけでは何もわからない。低気圧や高気圧が図—天気図—としてあらわれる。明日の天気がだいたいわかる。この天気図に入っている情報は、もともとラジオのアナウンサーが読み上げる数値データに全て含まれているものである。明日の天気を知るのはほとんど不可能である。紙の上に天気図という図にしてみてはじめてデータに隠されていた情報を引き出すことができる。[#h364e87b]データの可視化には可視化の目的や、対象データの種類の応じて様々な方法が開発されている。その中で、ここでは、スカラー場の可視化について簡単に説明する。[#za71e69f]xのある関数f(x)を直感的に理解するためには、グラフを書くのが一番である。xのx乗の関数、つまりf(x)=x^xはどんな形だろうか?下のグラフを見ればわかる。[#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を「挟めば」その辺を等高線が通る。辺上のどの位置を等高線が横切るかは、線形補間をすればよい。下の図はちょうど中点を通る例である。squaresと呼ばれる。等高線による可視化とは2次平面上に分布するスカラー場f(x,y)を等高線という曲線を見ることで理解する方法である。[#va1d6f14]等値面の「3次元版」も考えられる。この場合は3次元空間中の曲面の形状を見ることで3次元空間に分布する関数f(x,y,z)の分布を理解する。下の図の白い矢印で示している物体は等値面の例である。これは核融合実験装置(LHD)内部のプラズマの圧力分布の等値面である。Cubesアルゴリズムは、1985年に米国で特許となったが、20年を経過したので、現在では自由に利用できる。[#zd84152f]このような可視化アルゴリズムを実装したソフトウェアは既に各種開発されており、可視化ソフトウェアとかグラフソフトなどと呼ばれる。この演習では、そのなかで無償で利用出来るgnuplotを利用する。(なお、gnuplotはGNUのソフトウェアではない。/mnt/nfs/usr/bin/gnuplotというコマンド打てばgnuplotが立ち上がる。sin関数がプロットされるはずである。x**xとすれば見える。gnuplotの終了はコマンドプロンプトでquit(またはexit)と打つ。gnuplotではこのように解析的に定義された関数をターミナルに打ち込んでそのグラフを描く機能がある。それに加えて、ファイルに書き込まれた離散データを読み込み、それをグラフにする機能もある。この機能を使ってtemperature.f90で計算した正方形領域の中心点での温度が、ヤコビ法によって次第に収束して行く様子を可視化してみよう。heat1c.f90の出力(標準出力)をファイルtest.dataという名前に保存せよ。エディタで開くよりもheadコマンドやtailコマンドで見る方が早い。lpはlinespointsの略で、線(line)と点(point)を表示することを意味する。linespointsと書いてもよい。上の図の中で文字とグラフが重なって見にくいので、文字を消す。縦軸の表示範囲を調整をする。図全体のタイトルと、x軸、y軸の説明を入れる。[#fdc3f653]heat1c.f90で計算された最終的な平衡温度分布をgnuplotで見てみよう。正方形を真ん中で横に切るy=0.5の線上での温度のx分布をグラフにする。これからの演習でデータファイルを多数生成するので、データファイル出力先のディレクトリを用意しよう。../dataとしてdataディレクトリを作成する。別の場所、名前でも構わないが、以後のサンプルコードではそれに応じて適宜dataディレクトリの名前を変更すること。以下のプログラムをheat1_print_final_x_prof.f90という名前のソースコードにコピーし、コンパイル・実行せよ。heat1_print_final_x_profディレクトリ../dataにtemp.final_profile_xという名前のファイルが出来ているはずである。[#q7d70abf]../data/temp.final_profile_xのデータをgnuplotでグラフにする。このデータの中身をエディタまたはheadコマンド、tailコマンドで確認せよ。gnuplotではコマンドプロンプトに入力する内容をスクリプトファイルから読み込ませることが出来る。以下の内容をheat1_print_final_x_prof.gpという名前のファイルに保存せよ。heat1_print_final_x_prof.gpと打て。ファイルheat1_print_final_x_prof.gpを自由に変更してグラフがどう変わるか試せ。gnuplotのヘルプはコマンドプロンプトでhelpと打てば得られる。[#u4f89dc0]次にアニメーションによって収束の様子を確認しよう。そのためのデータ(連番つきファイル群)を書き出すためのプログラムheat_animation_x_prof.f90を以下に示す。実行すればdataディレクトリに連番ファイルが出力される。[#o7126c27]gnuplotのスクリプトを使えばアニメーションも簡単にできる。「gnuplotスクリプト自動生成プログラム」を作ったほうがよい。[#b2b4e3cd]以下のコードをheat1_animation_x_prof_gp_generator.f90というファイルに保存し、ソースコードの冒頭にあるUsageに従ってgnuplotによる可視化アニメーションを作れ。parallel__finalizeルーチン名を見れば意味は明白であろう。stopwatch__stt()に見られるように、関数名(サブルーチン名)は、「モジュール名」+「アンダースコア二つ」+「その関数の中身を示す言葉」で統一されている。その関数がどのモジュールで定義されているかがすぐに分かるようにするために工夫であった。[#k39ed7ce]以降では、2次元データの可視化を行う。これまでと同様gnuplotを利用する。[#a2717789]temperature.f90にある2次元データ出力ルーチンoutput_2d_profileは、正方形上(x,y平面上)に分布する温度をすべて書き出すものである。[#m11f0188]dataディレクトリ中の連番つきファイルtemp.2d.????の中身は以下のようになっているはず。data/temp.2d.0100のファイルに記された温度の分布をgnuplotの等高線で可視化してみよう。ファイル名やパラメータ等を自由に変更してその効果を試せ。時間発展(アニメーション)で等高線を見てみよう。例によってスクリプト生成プログラムを使う。利用方法は以下のコメント(Usage)に書かれている。温度を色で表現することも可能である。gnuplotのサンプルスクリプトは以下のとおり。時間発展のアニメーションを作ってみよう。2次元温度分布T(x,y)を高さ(z)で表すことも可能である。この時の描画は3次元的に行う必要がある。view)とも呼ばれる。gnuplotで描いてみよう。アニメーションを作ってみよう。gnuplotではviewというパラメータで視線の方向(「鳥」がどの方向から見下ろしているか)を指定することが出来る。このviewパラメータを変更したアニメーションを作ってみよう。0<=y<=1)が全体に一様に+4の強さで加熱されている場合を想定して計算していた。次に、下の図のように正方形を4つの部分にわけて市松模様状に加熱と冷却の分布を与えることを考える。冷却は負の加熱と考え、+4の加熱と-4の加熱の分布を図のように与える。答えの分布をsample_birdseyeview.gpで可視化すると次のような図が得られる。[#od9affb6]【提出は次週までにメールで(宛先:ymasada_at_harbor.kobe-u.ac.jp)】~熱伝導問題を解くコードを応用し、上に述べた市松模様状に+4と-4の加熱と冷却の分布を与えたときの熱平衡温度分布を求めるプログラムを作成せよ。市松模様状の加熱と冷却を与えるために、熱伝導問題を解くコードにどのような変更を加えたか説明せよ。市松模様状に加熱と冷却を与えた際の平衡解(最終状態解)温度の最大値(正の値)と最小値(負の値)を記せ。その平衡解の温度分布を可視化した図をつけること。gnuplotでなくても可。ファイルサイズは100kB以下にすること。[#veca5581](どれか一つ、一度だけ押してください。[#q4a8df41]「お名前」欄は空欄で可。