6.データ可視化 のバックアップ(No.2)
- バックアップ一覧
- 差分 を表示
- 現在との差分 を表示
- ソース を表示
- 6.データ可視化 へ行く。
- 1 (2010-07-07 (水) 10:20:45)
- 2 (2010-07-07 (水) 14:57:24)
- 3 (2010-07-07 (水) 19:16:50)
- 4 (2010-07-07 (水) 22:05:19)
- 5 (2010-07-07 (水) 23:07:31)
- 6 (2010-07-08 (木) 12:17:41)
- 7 (2010-07-14 (水) 17:46:39)
- 8 (2012-06-28 (木) 14:06:46)
- 9 (2013-07-02 (火) 18:23:15)
- 10 (2013-07-03 (水) 13:34:43)
- 11 (2013-07-03 (水) 16:03:59)
- 12 (2013-07-04 (木) 15:29:29)
- 13 (2013-07-11 (木) 12:56:31)
- 14 (2013-07-18 (木) 14:37:58)
- 15 (2014-05-08 (木) 14:15:12)
- 16 (2014-06-05 (木) 14:00:36)
現在（2014-06-11 (水) 22:33:33）作成中です。 既に書いている内容も大幅に変わる可能性が高いので注意。
神戸大学 大学院システム情報学研究科 計算科学専攻 陰山 聡
前回のレポート（並列化）の解説 †
問題設定の復習 †
#ref(): File not found: "problem_thermal.jpg" at page "6.データ可視化"
- 2次元正方形領域 [0,1]×[0,1] での熱伝導を考える
- 境界をすべて0℃に固定
- 領域全体に一定の熱を加える
- ⇒ このとき，十分な時間が経った後での温度分布はどうなるか？
非並列版コードの復習 †
熱伝導とは注目する一点が、「近所」の温度と等しくなろうという傾向である。 ⇒周囲の温度の平均値になろうとする。 熱源があれば、一定の割合で温度が上がろうとする。
ヤコビ法アルゴリズム。
#ref(): File not found: "jacobi_method.jpg" at page "6.データ可視化"
do j=1, m
do i=1, m
uij(n+1) = (ui-1,j(n) + ui+1,j(n) + ui,j-1(n) + ui,j+1(n) ) / 4 + fij
end do
end do
このヤコビ法に基づいたコード heat1.f90を前回の演習で見た。
heat1_wi_comments.f90 （詳細コメント付きバージョン） †
今回は、まずこのheat1.f90を少し丁寧に見る。
!
! heat1_wi_comments.f90
!
!-----------------------------------------------------------------------
! Time development form of the thermal diffusion equation is
! \partial T(x,y,t) / \partial t = \nabla^2 T(x,y,t) + heat_source.
! In the stationary state,
! \nabla^2 T(x,y) + heat_source = 0.
! The finite difference method with grid spacings dx and dy leads to
! (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.
! When 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.
! This suggests a relaxation method called Jacobi method adopted in
! this code.
!-----------------------------------------------------------------------
!
! heat1_wi_comments.f90
!
!-----------------------------------------------------------------------
! Time development form of the thermal diffusion equation is
! \partial T(x,y,t) / \partial t = \nabla^2 T(x,y,t) + heat_source.
! In the stationary state,
! \nabla^2 T(x,y) + heat_source = 0.
! The finite difference method with grid spacings dx and dy leads to
! (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.
! When 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.
! This suggests a relaxation method called Jacobi method adopted in
! this code.
!-----------------------------------------------------------------------
program heat1_wi_comments
implicit none
integer, parameter :: m=31 ! mesh size
integer :: nmax=20000 ! max loop counter
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 ! temperature field
real(DP), dimension(:,:), allocatable :: un ! work array
real(DP) :: heat=1.0_DP ! heat source term; uniform distribution.
real(DP) :: h ! grid size
! |<----- 1.0 ----->|
! j=m+1+-------u=0-------+ ---
! j=m | | ^
! . | | |
! . | uiform | |
! u=0 heat u=0 1.0
! j-direction . | source | |
!(y-direction) j=3 | | |
! ^ j=2 | | v
! | j=1 +-------u=0-------+ ---
! | i=0 1 2 ... i=m+1
! |
! +-------> i-direction (x-direction)
!
! when m = 7,
! |<-------------------- 1.0 -------------------->|
! | |
! | h h h h h h h h |
! |<--->|<--->|<--->|<--->|<--->|<--->|<--->|<--->|
! +-----+-----+-----+-----+-----+-----+-----+-----+
! i=0 1 2 3 4 5 6 7 8
!
if ( mod(m,2)==0 ) then
print *, 'm must be odd to have a grid on the center.'
stop
end if
allocate( u(0:m+1,0:m+1) ) ! memory allocation of 2-D array.
! you can access each element of u
! as u(i,j), where 0<=i<=m+1
! and 0<=j<=m+1.
allocate( un(m,m) ) ! another memory allocation. in
! this case un(i,j) with
! 1<=i<=m and 1<=j<=m.
h = 1.0_DP/(m+1) ! grid spacing.
u(:,:) = 0.0_DP ! initial temperature is zero all over the square.
do n = 1 , nmax ! relaxation loop
do j = 1 , m ! in the y-direction
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 ! same as the following do-loops:
u(1:m,1:m) = un(1:m,1:m) ! do j = 1 , m
! do i = 1 , m
! u(i,j) = un(i,j)
! end do
! end do
if ( mod(n,100)==0 ) print *, n, u(m/2+1,m/2+1) ! temperature at
end do ! the center
end program heat1_wi_comments
再びheat1.f90（コメントなしコード） †
次に詳しいコメントをはずしてもう一度もとのheat1.f90を見てみる。 （ただし、先週のコードから本質的でないところを少しだけ修正している。）
program heat1
implicit none
integer, parameter :: m=31, nmax=20000
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
if ( mod(m,2)==0 ) then
print *, 'm must be odd to have a grid on the center.'
stop
end if
allocate(u(0:m+1,0:m+1))
allocate(un(m,m))
h=1.0_DP/(m+1)
u(:,:) = 0.0_DP
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) print *, n, u(m/2+1,m/2+1)
end do
end program heat1
先週の課題（演習3-3） †
- heat1.f90 を MPI を用いて並列化せよ
解答例（heat2.f90） †
program heat2
use mpi
implicit none
integer, parameter :: m=31, nmax=20000
integer :: i,j,jstart,jend,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
integer :: nprocs,myrank,ierr,left,right
integer, dimension(MPI_STATUS_SIZE) :: istat
if ( mod(m,2)==0 ) then
print *, 'm must be odd to have a grid on the center.'
stop
end if
call mpi_init(ierr)
call mpi_comm_size(MPI_COMM_WORLD,nprocs,ierr)
call mpi_comm_rank(MPI_COMM_WORLD,myrank,ierr)
jstart=m*myrank/nprocs+1
jend=m*(myrank+1)/nprocs
allocate(u(0:m+1,jstart-1:jend+1))
allocate(un(m,jstart:jend))
!
! when m = 7, nprocs = 3
!
! |<-------------------- 1.0 -------------------->|
! | |
! | h | h h h h h h h |
! |<--->|<--->|<--->|<--->|<--->|<--->|<--->|<--->|
! | | | | | | | | |
! +-----+-----+-----+-----+-----+-----+-----+-----+
! j=0 1 2 3 4 5 6 7 8
! +-----+-----+-----+-----+-----+-----+-----+-----+
! | jstart | | | | | |
! | | jend | | jstart | |
! o-----rank=0------o | | jend |
! | | | | | |
! | jstart jend | | |
! | | | | | |
! o------- rank=1---------o | |
! | | | |
! o------rank=3-----o
!
h = 1.0_DP/(m+1)
u(:,:) = 0.0_DP
left=myrank-1
if (myrank==0) left=nprocs-1
right=myrank+1
if (myrank==nprocs-1) right=0
do n=1, nmax
call mpi_sendrecv(u(1,jend),m,MPI_DOUBLE_PRECISION,right,100, &
& u(1,jstart-1),m,MPI_DOUBLE_PRECISION,left,100, &
& MPI_COMM_WORLD,istat,ierr)
call mpi_sendrecv(u(1,jstart),m,MPI_DOUBLE_PRECISION,left,100, &
& u(1,jend+1),m,MPI_DOUBLE_PRECISION,right,100, &
& MPI_COMM_WORLD,istat,ierr)
if (myrank==0) u(1:m,0)=0.0_DP
if (myrank==nprocs-1) u(1:m,m+1)=0.0_DP
do j=jstart, jend
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,jstart:jend)=un(1:m,jstart:jend)
if (jstart<=m/2+1 .and. jend>=m/2+1) then
if (mod(n,100)==0) print *, n, u(m/2+1,m/2+1)
end if
end do
call mpi_finalize(ierr)
end program heat2
コード解説 （コメント付きソースコード heat2_wi_comments.f90） †
上のheat2.f90を解説するために詳しいコメントをつけたソースコードheat2_wi_comments.f90を以下に示す。
program heat2_wi_comments
use mpi
implicit none
integer, parameter :: m=31 ! mesh size
integer, parameter :: nmax=20000 ! max loop counter
integer :: i,j,jstart,jend,n
integer, parameter :: SP = kind(1.0)
integer, parameter :: DP = selected_real_kind(2*precision(1.0_SP))
real(DP), dimension(:,:), allocatable :: u ! temperature field
real(DP), dimension(:,:), allocatable :: un ! work array
real(DP) :: heat=1.0_DP ! heat source term. uniform distribution.
real(DP) :: h ! grid spacing; dx=dy=h
integer :: nprocs ! number of MPI processes
integer :: myrank ! my rank number
integer :: left, right ! nearest neighbor processes
integer, dimension(MPI_STATUS_SIZE) :: istat ! used for MPI
integer :: ierr ! used for mpi routines
if ( mod(m,2)==0 ) then ! to print out the centeral temperature.
print *, 'm must be odd to have a grid on the center.'
stop
end if
call mpi_init(ierr)
call mpi_comm_size(MPI_COMM_WORLD,nprocs,ierr)
call mpi_comm_rank(MPI_COMM_WORLD,myrank,ierr)
h = 1.0_DP/(m+1) ! grid spacings. see the comment fig below.
jstart = m*myrank/nprocs+1 ! you are in charge of this point
jend = m*(myrank+1)/nprocs ! ...to this points. see the fig below.
!
! when m = 7, nprocs = 3
!
! |<-------------------- 1.0 -------------------->|
! | |
! | h | h h h h h h h |
! |<--->|<--->|<--->|<--->|<--->|<--->|<--->|<--->|
! | | | | | | | | |
! +-----+-----+-----+-----+-----+-----+-----+-----+
! j=0 1 2 3 4 5 6 7 8
! +-----+-----+-----+-----+-----+-----+-----+-----+
! | jstart | | | | | |
! | | jend | | jstart | |
! o-----rank=0------o | | jend |
! | | | | | |
! | jstart jend | | |
! | | | | | |
! o------- rank=1---------o | |
! | | | |
! o------rank=3-----o
!
allocate( u(0:m+1,jstart-1:jend+1)) ! 2-D array for temperature.
allocate(un(m,jstart:jend)) ! 2-D work array.
u(:,:) = 0.0_DP
left = myrank-1 ! left neighbor.
if (myrank==0) left = nprocs-1 ! a dummy periodic distribution
right = myrank+1 ! right neighbor. ! for code simplicity. this has
if (myrank==nprocs-1) right = 0 ! actually no effect.
do n = 1 , nmax ! main loop for the relaxation, to equilibrium.
call mpi_sendrecv(u(1,jend), m,MPI_DOUBLE_PRECISION,right,100, &
u(1,jstart-1),m,MPI_DOUBLE_PRECISION,left,100, &
MPI_COMM_WORLD,istat,ierr)
call mpi_sendrecv(u(1,jstart),m,MPI_DOUBLE_PRECISION,left, 100, &
u(1,jend+1),m,MPI_DOUBLE_PRECISION,right,100, &
MPI_COMM_WORLD,istat,ierr)
if (myrank==0) u(1:m,0 ) = 0.0_DP ! to keep the boundary
if (myrank==nprocs-1) u(1:m,m+1) = 0.0_DP ! condition (temp=0).
do j = jstart , jend ! Jacobi method.
do i = 1 , m ! no need to calculate the boundary i=0 and m+1.
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,jstart:jend)=un(1:m,jstart:jend) ! this is actually doubled do-loops.
if (jstart<=m/2+1 .and. jend>=m/2+1) then ! print out the temperature
if (mod(n,100)==0) print *, n, u(m/2+1,m/2+1) ! at the central grid point.
end if
end do
deallocate(un,u) ! or you can call deallocate for two times for un and u.
call mpi_finalize(ierr)
end program heat2_wi_comments
可視化の必要性 †
可視化とは †
ここではスカラーデータの可視化のみ。
1次元グラフ †
2次元等高線アルゴリズム †
3次元等値面アルゴリズム †
コードのリファクタリング †
今後の演習で行うコードの改訂
- 1次元可視化機能の追加（中心点の温度の収束の様子を示す）
- 2次元可視化機能の追加（温度場全体の収束の様子をアニメーションで示す）
- 2次元領域分割による並列化
- 問題の3次元化
gnuplot入門 †
1次元可視化（中心点の温度の収束の様子を示す） †
【演習】 †
2次元可視化（温度場全体の収束の様子のアニメーション） †
【演習】 †
as of 2024-04-27 (土) 08:11:49 (5141)