現在(2015-11-07 (土) 15:06:03)作成中です。 既に書いている内容も大幅に変わる可能性が高いので注意。
神戸大学 大学院システム情報学研究科 計算科学専攻 陰山 聡
【目次】
pause -1と書く必要があった。
前回の宿題は、「thermal_diffusion.f90を理解しておくこと」であった。 ここではポイントだけ説明する。
モジュール名 | 内容 |
constants | 定数 |
parallel | 並列化 |
temperature | 温度場 |
データ出力ルーチンを除くとコードの構造は以下のとおり:
call parallel__initialize call temperature__initialize do loop = 1 , LOOP_MAX call temperature__update end do call temperature__finalize call parallel__finalize
ルーチン名を見れば意味は明白であろう。
call parallel__initialize call temperature__initialize call temperature__output_1d_profile call temperature__output_2d_profile
に見られるように、関数名(サブルーチン名)は、「モジュール名」+「アンダースコア二つ」+「その関数の中身を示す言葉」で統一されている。 その関数がどのモジュールで定義されているかがすぐに分かるようにするために工夫である。
今考えている正方形領域は 一様かつ一定に加熱されている。 その加熱率 heat_source を(ある規格化の下で)4とすると、 定常状態における温度T(x,y)は以下のポアッソン方程式を満たす。
つまり
である。
連続場としての温度場を 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を繰り返し計算する。
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という定数を足すだけで済むようにしている。
thermal_diffusion.f90をコンパイルし、計算機「scalar」で実行せよ。
今回は2次元データの可視化を行う。 前回と同様gnuplotを利用する。
gnuplotで2次元データをプロットする場合、そのデータのフォーマットは
x00 y00 関数値 x01 y00 関数値 x02 y00 関数値 ・ ・ ・ x09 y00 関数値 (空行) x00 y01 関数値 x01 y01 関数値 x02 y01 関数値 ・ ・ ・ x09 y01 関数値 (空行) ・ ・ ・ x09 y09 関数値
である。
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
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 ・ ・ ・
静止画による説明
静止画による説明
静止画による説明
ソースコードの変更。 thermal_diffusion2.f90
thermal_diffusion02.f90で求めたチェッカーボード型熱源分布の平衡解(最終状態解)を可視化せよ。 温度の最大値と最小値も記せ。
今回の演習内容はどうでしたか?(どれか一つ、一度だけ押してください。)
「お名前」欄は空欄で可。
最新の10件を表示しています。 コメントページを参照
as of 2024-03-19 (火) 19:14:58 (6100)