PukiWiki - 6.データ可視化のバックアップ(No.3)

現在（2014-06-11 (水) 22:33:33）作成中です。既に書いている内容も大幅に変わる可能性が高いので注意。

神戸大学大学院システム情報学研究科計算科学専攻陰山聡

前回のレポート（並列化）の解説
サンプルジョブスクリプト
【演習】
- 結果例
可視化の必要性
可視化の例
熱伝導問題コードのリファクタリング
thermal_diffusion.f90コード
gnuplot入門
1次元可視化（中心点の温度の収束の様子を示す）
【演習】
2次元可視化（温度場全体の収束の様子のアニメーション）
【演習】
授業アンケート

前回のレポート（並列化）の解説 †

↑

問題設定の復習 †

#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)

  h = 1.0_DP/(m+1)

  jstart = m*myrank/nprocs+1
  jend   = m*(myrank+1)/nprocs

  allocate( u(0:m+1,jstart-1:jend+1))
  allocate(un(    m,  jstart:jend  ))

  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=2-----o
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
! Another example, when m=17, nprocs=6
!
!  j=0  1  2  3  4  5  6  7  8  9  10 11 12 13 14 15 16 17 18
!    +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
!    |        |        |     |           |        |        |
!    o==rank0=o        |     o===rank3===o        |        |
!          |           |        |     |           |        |
!          o===rank1===o        |     o===rank4===o        |
!                   |           |              |           |
!                   o===rank2===o              o===rank5===o
! each process      |           |
! has 2D array  u[i,jstart-1:jend+1] 
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  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.
    !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !
    !                  jend                          jstart      
    !   <(left)>        |                              |    <(right)>
    !      - --o--o--o--o--o                        o--o--o--o--o-- -
    !              send |  ^ recv                   ^  |
    !                   |  |                        |  |
    !                 [2]  [3]                    [1]  [4]
    !                   |  |                        |  |
    !              recv v  | send                   |  v
    !                   o--o--o--o-- - - --o--o--o--o--o
    !                      |                        |
    !                    jstart                    jend
    !                               <(myself)>
    !                     
    !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    call mpi_sendrecv(u(1,jend),    m,MPI_DOUBLE_PRECISION,right,100,  & ! [1]
                      u(1,jstart-1),m,MPI_DOUBLE_PRECISION,left,100,   & ! [2]
                      MPI_COMM_WORLD,istat,ierr)
    call mpi_sendrecv(u(1,jstart),m,MPI_DOUBLE_PRECISION,left, 100,    & ! [3]
                      u(1,jend+1),m,MPI_DOUBLE_PRECISION,right,100,    & ! [4]
                      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

↑

MPI通信の確認 (heat2_sendrecv_check.f90) †

前回の演習で説明されたMPI_SENDRECVの挙動を改めて確認しよう。そのためのソースコード（heat2_sendrecv_check.f90）を以下に用意した。これは上記のheat2.f90からヤコビ法による熱伝導問題の解法部分を除き、 MPIの通信に関係する部分だけを抜き出したものである。

!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
! 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=2 ====o
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
! Another example, when m=17, nprocs=6
!
!  j=0  1  2  3  4  5  6  7  8  9  10 11 12 13 14 15 16 17 18
!    +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
!    |        |        |     |           |        |        |
!    o==rank0=o        |     o===rank3===o        |        |
!          |           |        |     |           |        |
!          o===rank1===o        |     o===rank4===o        |
!                   |           |              |           |
!                   o===rank2===o              o===rank5===o
! each process      |           |
! has 2D array  u[i,jstart-1:jend+1] 
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  
program heat2_sendrecv_check
  use mpi
  implicit none
  integer, parameter :: m=17
  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
  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  ))

  u(:,jstart-1:jend+1) = myrank

  call print('before',myrank, jstart, jend, u(m/2+1,:))

  left = myrank-1
  if ( myrank==0 )        left  = nprocs-1
  right = myrank+1
  if ( myrank==nprocs-1 ) right = 0

  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !
  !                  jend                          jstart      
  !   <(left)>        |                              |    <(right)>
  !      - --o--o--o--o--o                        o--o--o--o--o-- -
  !              send |  ^ recv                   ^  |
  !                   |  |                        |  |
  !                 [2]  [3]                    [1]  [4]
  !                   |  |                        |  |        
  !                   |  |                        |  |
  !              recv v  | send                   |  v
  !                   o--o--o--o-- - - --o--o--o--o--o
  !                      |                        |
  !                    jstart                    jend
  !                               <(myself)>
  !                     
  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  call mpi_sendrecv(u(1,jend),    m,MPI_DOUBLE_PRECISION,right,100,   & ! [1]
                    u(1,jstart-1),m,MPI_DOUBLE_PRECISION, left,100,   & ! [2]
                    MPI_COMM_WORLD,istat,ierr)
  call mpi_sendrecv(u(1,jstart),m,MPI_DOUBLE_PRECISION, left,100,     & ! [3]
                    u(1,jend+1),m,MPI_DOUBLE_PRECISION,right,100,     & ! [4]
                    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

  call print(' after', myrank, jstart, jend, u(m/2+1,:))

  deallocate(u,un)
  call mpi_finalize(ierr)

contains

  subroutine print(message,myrank,jstart,jend,u)
    character(len=*), intent(in)                     :: message
    integer,          intent(in)                     :: myrank
    integer,          intent(in)                     :: jstart
    integer,          intent(in)                     :: jend
    real(DP), dimension(jstart-1:jend+1), intent(in) :: u

    integer, dimension(0:m+1) :: iu          ! work array of integer
    integer                   :: rank, ierr
    character(len=100)        :: line        ! printed line on stdout
    character(len=20)         :: format      ! print format
    
    iu(:) = -1    ! initialization of the work integer array.
    iu(jstart-1:jend+1) = u(jstart-1:jend+1) ! copy the data (double to int).

    write(format,'(i3"(i3)")') m+2     ! you are not expected to 
    write(line,format) (iu(i),i=0,m+1) ! understand this trick.

    do rank = 0 , nprocs-1                   ! to avoid random printing
       call mpi_barrier(MPI_COMM_WORLD,ierr) ! by each process in random order,
       if ( rank==myrank ) then              ! we here put the loop and barrier.
          print *, message//':'//trim(line)
       end if
    end do
  end subroutine print

end program heat2_sendrecv_check

このコード（heat2_sendrecv_check.f90）をscalarに投入するためのジョブスクリプトを以下のようにすると、このコードのコメントにあるとおり、グリッドサイズmが17で、プロセス数が6のジョブを投入できる。

↑

サンプルジョブスクリプト †

#PBS -l cputim_job=00:05:00
#PBS -l memsz_job=2gb
#PBS -l cpunum_job=1
#PBS -T vltmpi
#PBS -b 6
#PBS -q PCL-A

cd /home/users/yourID/YourDirectory
mpirun_rsh -np 6 ${NQSII_MPIOPTS} ./a.out

これをたとえば heat2_sendrecv_check.jsという名前のファイルに保存して

qsub heat2_sendrecv_check.js

とすればジョブが投入できる。

↑

【演習】 †

上のheat2_sendrecv_check.f90をマシン「scalar」で計算せよ。

↑

結果例 †

%NQSII(INFO): ------- Output from job:0000 -------
[0]  before:  0  0  0  0 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
[0]   after:  5  0  0  1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
%NQSII(INFO): ------- Output from job:0001 -------
[1]  before: -1 -1  1  1  1  1  1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
[1]   after: -1 -1  0  1  1  1  2 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
%NQSII(INFO): ------- Output from job:0002 -------
[2]  before: -1 -1 -1 -1 -1  2  2  2  2  2 -1 -1 -1 -1 -1 -1 -1 -1 -1
[2]   after: -1 -1 -1 -1 -1  1  2  2  2  3 -1 -1 -1 -1 -1 -1 -1 -1 -1
%NQSII(INFO): ------- Output from job:0003 -------
[3]  before: -1 -1 -1 -1 -1 -1 -1 -1  3  3  3  3  3 -1 -1 -1 -1 -1 -1
[3]   after: -1 -1 -1 -1 -1 -1 -1 -1  2  3  3  3  4 -1 -1 -1 -1 -1 -1
%NQSII(INFO): ------- Output from job:0004 -------
[4]  before: -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1  4  4  4  4  4 -1 -1 -1
[4]   after: -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1  3  4  4  4  5 -1 -1 -1
%NQSII(INFO): ------- Output from job:0005 -------
[5]  before: -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1  5  5  5  5  5
[5]   after: -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1  4  5  5  5  0

↑

可視化の必要性 †

上の例では、MPIプロセスがもつ配列uの値をテキストデータとして書き出すことで、プロセス間のMPI_SENDRECVによる通信の様子が分かりやすく表示された。このように人間の視覚を使って情報を引き出す操作を一般に可視化という。

科学技術計算においては、可視化が様々な局面で重要な役割を果たす。特に、計算結果として出力される数値データに潜む情報を引き出すための可視化―これをデータ可視化という―は、不可欠と言ってよい。

データ可視化の身近な例は天気図であろう。

NKHラジオの「気象通報」という番組を知っているだろうか？各地の風向、風力、気圧などの数字をアナウンサーが次々に読み上げていくラジオ番組である。これをただ耳で聞いているだけでは何もわからない。しかしペンと地図を手にし、番組を聞きながら読み上げられる各地点に気圧などを書き込み（専用の地図を売っている）、最後に同じ気圧になりそうな点を曲線（等高線）でつなげていくと、低気圧や高気圧が図―天気図―としてあらわれる。知識のある人が天気図を見れば、明日の天気がだいたいわかる。

この天気図に入っている情報は、もともとラジオのアナウンサーが読み上げる数値データに全て含まれているものである。だが、そのような数値データの羅列をいくら注意深く聞いても、あるいはそれを数表にして紙の上に書き、それをじっと睨んでみても、明日の天気を知るのはほとんど不可能である。紙の上に天気図という図にしてみてはじめてデータに隠されていた情報を引き出すことができる。これが可視化である。

↑

可視化の例 †

データの可視化には可視化の目的や、対象データの種類の応じて様々な方法が開発されている。その中で、ここでは、スカラー場の可視化について簡単に説明する。

↑

1次元グラフ †

xのある関数f(x)を直感的に理解するためには、グラフを書くのが一番である。例えばxのx乗の関数、つまりf(x)=x^xはどんな形だろうか？下のグラフを見ればわかる。

#ref(): File not found: "graph_1d.png" at page "6.データ可視化"

↑

2次元等高線 †

xとyの関数、つまり2次元の関数f(x,y)の形を理解するには、等高線を描くのがよい。地表面での大気の圧力pの分布p(x,y)の等高線は天気図でおなじみである。

#ref(): File not found: "graph_2d.png" at page "6.データ可視化"

上の図では等高線の高さを色で表示してさらに分かりやすくしている。

ここで等高線の描画アルゴリズムを紹介しよう。計算格子点上に定義されたデータから一本の等高線を描くには、下の図のように短い線分をつなげていけばよい。

#ref(): File not found: "contour_line01.png" at page "6.データ可視化"

一つの線分は4つの計算格子点で定義された長方形領域（セル）の中で直線を描く。例えばf(x,y)=1.0の値の等高線を描く場合を考えよう。あるセルの4つの頂点におけるfの値が全て1.0よりも大きいか、あるいは1.0未満であれば、f=1.0の等高線はこのセルを通らないのは明らかである。

セルを周囲の4つの辺、それぞれの両端の頂点でのfの値が1.0を「挟めば」その辺を等高線が通る。辺上のどの位置を等高線が横切るかは、線形補間をすればよい。下の図はちょうど中点を通る例である。

#ref(): File not found: "contour_line02.png" at page "6.データ可視化"

このアルゴリズムは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）内部のプラズマの圧力分布の等値面である。

#ref(): File not found: "isosurface_sample.png" at page "6.データ可視化"

このような等値面を描くためのアルゴリズムとしてMarching Cubesがある。

#ref(): File not found: "marching_cubes.png" at page "6.データ可視化"

このMarching Cubesアルゴリズムは、 1985年に米国で特許となったが、20年を経過したので、現在では自由に利用できる。

↑

熱伝導問題コードのリファクタリング †

可視化アルゴリズムの紹介はこれくらいにとどめて、これから前回の講義で学んだ2次元正方形領域での熱伝導問題を例にとり、簡単な可視化を実践的に演習してみよう。目標は、heat1.f90コードと、heat2.f90に

1次元可視化機能の追加（中心点の温度が収束する様子をグラフで示す）
2次元可視化機能の追加（温度場全体が収束していく様子をアニメーションで示す）

という可視化機能を追加することである。ただし、グラフや等高線などの可視化にはここではgnuplotというフリーウェアを利用する。

また、これまで学んだMPIによる並列化は、 1方向のみの分割による並列化（1次元領域分割による並列化）だけであったが、次回以降の演習で、もう少し実際的な並列化問題を扱う。具体的には、

2次元熱伝導問題の2次元領域分割による並列化
熱伝導問題の3次元化

を行う予定である。

したがって、今後扱うコードは少しづつ複雑化していくので、ここでまずはコードのリファクタリングをしておこう。

↑

thermal_diffusion.f90コード †

!
! thermal_diffusion.f90
!
module constants
  implicit none
  integer, parameter :: SP = kind(1.0)
  integer, parameter :: DP = selected_real_kind(2*precision(1.0_SP))
  integer, parameter :: MESH_SIZE = 61
  integer, parameter :: LOOP_MAX = 10000
end module constants

module parallel
  use constants
  use mpi
  implicit none
  private
  public :: parallel__initialize,       &
            parallel__finalize,         &
            parallel__sendrecv,         &
            parallel__set_prof_1d,      &
            parallel__set_prof_2d,      &
            parallel__tellme
  
  integer :: nprocs
  integer :: myrank, left, right
  integer :: jstart, jend

contains

  subroutine parallel__finalize
    integer :: ierr
    call mpi_finalize(ierr)
  end subroutine parallel__finalize
  
  subroutine parallel__initialize
    integer :: ierr
    call mpi_init(ierr)
    call mpi_comm_size(MPI_COMM_WORLD, nprocs, ierr)
    call mpi_comm_rank(MPI_COMM_WORLD, myrank, ierr)
    left  = myrank - 1
    right = myrank + 1
    if ( myrank==0 )        left  = nprocs-1
    if ( myrank==nprocs-1 ) right = 0
    jstart = MESH_SIZE * myrank / nprocs + 1
    jend   = MESH_SIZE * (myrank+1) / nprocs
  end subroutine parallel__initialize

  subroutine parallel__sendrecv(send_vector,rank_target,        &
                                recv_vector,rank_source)
    real(DP), dimension(:), intent(in)  :: send_vector
    integer,                intent(in)  :: rank_target
    real(DP), dimension(:), intent(out) :: recv_vector
    integer,                intent(in)  :: rank_source

    integer, dimension(MPI_STATUS_SIZE) :: istat
    integer :: tag_dummy = 100
    integer :: n, ierr
    n = size(send_vector,dim=1)
    call mpi_sendrecv(send_vector,              &
                      n,                        &
                      MPI_DOUBLE_PRECISION,     & 
                      rank_target,              &
                      tag_dummy,                &
                      recv_vector,              &
                      n,                        &
                      MPI_DOUBLE_PRECISION,     &
                      rank_source,              &
                      tag_dummy,                &
                      MPI_COMM_WORLD,           &
                      istat,                    &
                      ierr)
  end subroutine parallel__sendrecv

  function parallel__set_prof_1d(jstt_,jend_,myprof) result(prof_1d)
    integer,  intent(in)                         :: jstt_
    integer,  intent(in)                         :: jend_
    real(DP), dimension(jstt_:jend_), intent(in) :: myprof
    real(DP), dimension(0:MESH_SIZE+1)           :: prof_1d

    real(DP), dimension(0:MESH_SIZE+1) :: work
    integer                            :: ierr

    work(:) = 0.0_DP
    work(jstt_:jend_) = myprof(jstt_:jend_)
    call mpi_allreduce(work,                    &  ! source
                       prof_1d,                 &  ! target
                       MESH_SIZE+2,             &
                       MPI_DOUBLE_PRECISION,    &
                       MPI_SUM,                 &
                       MPI_COMM_WORLD,          &
                       ierr)
  end function parallel__set_prof_1d

  function parallel__set_prof_2d(jstt_,jend_,myprof) result(prof_2d)
    integer,  intent(in)                             :: jstt_
    integer,  intent(in)                             :: jend_
    real(DP), dimension(0:MESH_SIZE+1,jstt_:jend_),                    &
                                          intent(in) :: myprof
    real(DP), dimension(0:MESH_SIZE+1,0:MESH_SIZE+1) :: prof_2d

    real(DP), dimension(0:MESH_SIZE+1,0:MESH_SIZE+1) :: work
    integer :: ierr
    integer :: meshsize_p1_sq = (MESH_SIZE+1)**2

    work(:,:) = 0.0_DP
    work(:,jstt_:jend_) = myprof(:,jstt_:jend_)
    call mpi_allreduce(work(1,1),               &  ! source
                       prof_2d(1,1),            &  ! target
                       meshsize_p1_sq,          &
                       MPI_DOUBLE_PRECISION,    &
                       MPI_SUM,                 &
                       MPI_COMM_WORLD,          &
                       ierr)
  end function parallel__set_prof_2d
  
  function parallel__tellme(which) result(val)
    character(len=*), intent(in) :: which
    integer                      :: val

    select case (which)
      case  ('jend')
        val = jend
      case  ('jstart')
        val = jstart
      case  ('left')
        val = left
      case  ('right')
        val = right
      case  ('myrank')
        val = myrank
      case  ('nprocs')
        val = nprocs
    end select
  end function parallel__tellme

end module parallel


module temperature
  use constants
  use parallel
  implicit none
  private 
  public :: temperature__initialize,            &
            temperature__finalize,              &
            temperature__output_1d_profile,     &
            temperature__output_2d_profile,     &
            temperature__update

  real(DP), allocatable, dimension(:,:) :: temp
  real(DP), allocatable, dimension(:,:) :: work
  real(DP), parameter :: SIDE = 1.0_DP
  real(DP) :: h = SIDE / MESH_SIZE
  real(DP) :: heat
  integer  :: jstart, jend
  integer  :: myrank, right, left, nprocs
  
  !   
  !               |<----------- SIDE ---------->|
  !               |                             |
  !               |<-h->|                       |
  !               0     1     2    ...       MESH_SIZE+1
  !        jend+1 +-----+-----+-----+-----+-----+
  !               |     |     |     |     |     |
  !               |     |     |     |     |     |
  !          jend +-----+-----+-----+-----+-----+  
  !               |     |     |     |     |     |  
  !            .                                   
  !            .      temp(0:MESH_SIZE,0:jend+1)   
  !            .                                   
  !               |     |     |     |     |     |  
  !        jstart +-----+-----+-----+-----+-----+  
  !               |     |     |     |     |     |    
  !               |     |     |     |     |     |
  !      jstart-1 +-----+-----+-----+-----+-----+
  ! 

contains

!----------------
!--- private  ---
!----------------

  subroutine set_jstart_and_jend(jstart_, jend_)
    integer, intent(out) :: jstart_, jend_
    if ( myrank==0 ) then
       jstart_ = 0
       jend_   = jend
    else if ( myrank==nprocs-1 ) then
       jstart_ = jstart
       jend_   = MESH_SIZE+1
    else
       jstart_ = jstart
       jend_   = jend
    end if
  end subroutine set_jstart_and_jend

!--------------
!--- public ---
!--------------

  subroutine temperature__initialize
    real(DP) :: heat_source = 4.0
    jstart = parallel__tellme('jstart')
    jend   = parallel__tellme('jend')
    myrank = parallel__tellme('myrank')
    right  = parallel__tellme('right')
    left   = parallel__tellme('left')
    nprocs = parallel__tellme('nprocs')
    allocate(temp(0:MESH_SIZE+1,jstart-1:jend+1))
    allocate(work(    MESH_SIZE,  jstart:jend)  )
    heat = (heat_source/4) * h * h
    temp(:,:) = 0.0_DP       ! initial condition
  end subroutine temperature__initialize

  subroutine temperature__finalize
    deallocate(work,temp)
  end subroutine temperature__finalize

  subroutine temperature__output_1d_profile
    real(DP), dimension(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.i=middle."
    integer, parameter          :: icut = MESH_SIZE/2+1 ! MESH_SIZE is odd
    !   nprocs = 2  
    !   MESH_SIZE = 9
    ! 
    !    0     1     2     3     4     5     6     7     8     9    10
    !    +-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+
    !    |<=====rank0=====>|     |<==rank1==>|     |<=====rank1=====>|
    !  jstart_ |         jend_   |           |     |           |     |
    !    |   jstart      jend  jstart_     jend_   |           |     |
    !    |     |           |   jstart      jend  jstart_       |   jend_  
    !    |     |           |     |           |   jstart      jend    |
    !    |     |           |     |           |     |           |     |
    !    +-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+
    !
    integer :: j
    call set_jstart_and_jend(jstart_,jend_)
    write(serial_num,'(i4.4)') counter
    prof(:) = parallel__set_prof_1d(jstart_,jend_,temp(icut,jstart_:jend_))
    if ( myrank==0 ) then
       open(10,file=base//serial_num)
       do j = 0 , MESH_SIZE+1
          write(10,*) j, prof(j)
       end do
       close(10)
    end if
    counter = counter + 1
  end subroutine temperature__output_1d_profile

  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."
    !   nprocs = 2  
    !   MESH_SIZE = 9
    !    0     1     2     3     4     5     6     7     8     9    10
    !    +-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+
    !    |<=====rank0=====>|     |<==rank1==>|     |<=====rank1=====>|
    !  jstart_ |         jend_   |           |     |           |     |
    !    |   jstart      jend  jstart_     jend_   |           |     |
    !    |     |           |   jstart      jend  jstart_       |   jend_  
    !    |     |           |     |           |   jstart      jend    |
    !    |     |           |     |           |     |           |     |
    !    +-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+
    !
    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

  subroutine temperature__update
    integer :: i, j
    call parallel__sendrecv(temp(1:MESH_SIZE,jend),     right,   &
                            temp(1:MESH_SIZE,jstart-1), left  )
    call parallel__sendrecv(temp(1:MESH_SIZE,jstart),   left,    &
                            temp(1:MESH_SIZE,jend+1),   right )
    if ( myrank==0 )        temp(1:MESH_SIZE,0          ) = 0.0_DP
    if ( myrank==nprocs-1 ) temp(1:MESH_SIZE,MESH_SIZE+1) = 0.0_DP
    !
    ! 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 by the following 
    ! update procedure called Jacobi method;
    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
    temp(1:MESH_SIZE,jstart:jend) = work(1:MESH_SIZE,jstart:jend)
  end subroutine temperature__update

end module temperature


program thermal_diffusion
  use constants
  use parallel
  use temperature
  implicit none

  integer :: loop

  call check
  call parallel__initialize
  call temperature__initialize
  call temperature__output_1d_profile
  call temperature__output_2d_profile
  do loop = 1 , LOOP_MAX
     call temperature__update
     if ( mod(loop,100)==0 ) call temperature__output_1d_profile
     if ( mod(loop,100)==0 ) call temperature__output_2d_profile
  end do
  call temperature__finalize
  call parallel__finalize

contains

  subroutine check
    if ( mod(MESH_SIZE,2)==0 ) then
       print *, 'MESH_SIZE must be odd to have a grid on the center.'
       stop
    end if
  end subroutine check
end program thermal_diffusion

↑

gnuplot入門 †

↑

1次元可視化（中心点の温度の収束の様子を示す） †

↑

【演習】 †

↑

2次元可視化（温度場全体の収束の様子のアニメーション） †

↑

【演習】 †

↑

授業アンケート †

as of 2024-04-19 (金) 07:42:24 (5141)

6.データ可視化のバックアップ(No.3)

前回のレポート（並列化）の解説 †

問題設定の復習 †

非並列版コードの復習 †

heat1_wi_comments.f90 （詳細コメント付きバージョン） †

再びheat1.f90（コメントなしコード） †

先週の課題（演習3-3） †

解答例（heat2.f90） †

コード解説（コメント付きソースコード heat2_wi_comments.f90） †

MPI通信の確認 (heat2_sendrecv_check.f90) †

サンプルジョブスクリプト †

【演習】 †

結果例 †

可視化の必要性 †

可視化の例 †

1次元グラフ †

2次元等高線 †

3次元等値面 †

熱伝導問題コードのリファクタリング †

thermal_diffusion.f90コード †

gnuplot入門 †

1次元可視化（中心点の温度の収束の様子を示す） †

【演習】 †

2次元可視化（温度場全体の収束の様子のアニメーション） †

【演習】 †

授業アンケート †

演習について

計算科学演習A1

計算科学演習A2

πコンピュータ利用ガイド

このサイトについて

お問い合わせ

サイト内検索

6.データ可視化 のバックアップ(No.3)

前回のレポート（並列化）の解説 †

問題設定の復習 †

非並列版コードの復習 †

heat1_wi_comments.f90 （詳細コメント付きバージョン） †

再びheat1.f90（コメントなしコード） †

先週の課題（演習3-3） †

解答例（heat2.f90） †

コード解説 （コメント付きソースコード heat2_wi_comments.f90） †

MPI通信の確認 (heat2_sendrecv_check.f90) †

サンプルジョブスクリプト †

【演習】 †

結果例 †

可視化の必要性 †

可視化の例 †

1次元グラフ †

2次元等高線 †

3次元等値面 †

熱伝導問題コードのリファクタリング †

thermal_diffusion.f90コード †

gnuplot入門 †

1次元可視化（中心点の温度の収束の様子を示す） †

【演習】 †

2次元可視化（温度場全体の収束の様子のアニメーション） †

【演習】 †

授業アンケート †

演習について

計算科学演習A1

計算科学演習A2

πコンピュータ利用ガイド

このサイトについて

お問い合わせ

サイト内検索

6.データ可視化のバックアップ(No.3)

コード解説（コメント付きソースコード heat2_wi_comments.f90） †