メニーコアと格闘ブログ:ポアソン方程式の並列化 その3

出典: トータル・ディスクロージャ・サイト(事実をありのままに)

[前回まで]

1次元N個の格子点の計算を行うにあたり、1からNまで順番に計算すると依存関係が発生するため、複数の格子を同時に計算できないことを示しました。


これを踏まえて、偶数番目、奇数番目の格子点を交互に計算するよう、“計算の順番を操作する”ことで、計算に並列性を持たせることができることを示しました。

実はこの手法は“ゼブラ法”と呼ばれる手法で、以下に詳しいです。

並列プログラミング入門 MPI版

理化学研究所 情報基盤センター

2007牛:2月7日版

青山幸也

http://accc.riken.jp/HPC/training/mpi/mpi_all_2007-02-07.pdf (プログラム例(ゼブラ法):P5-59)


今回は、この並列性を使って、ポアソン方程式を並列に解く実装例を示します。

元々のソースコードは、ポアソン方程式をガウスザイデル法を使って、反復計算により解を特定する方法です。

[1プロセッサで計算する元のソースコード]

   do kk = 1,km
       g1 = 0.

       !*** boundary gondi'iioh‘for pressure
       !****y-z
       do k = 1,k21
       do j = 1,j21
           p(1,j,k) = p(2,j,k)
           p(i20,j,k) = p(i19,j,k)
       enddo
       enddo
       !***x-z
       do k = 1,k21
       do i = 1,i21
           p(i,1,k) = p(i,2,k)
           p(i,j20,k) = p(i,j19,k)
       enddo
       enddo
       !*** x-y
       do j = 1,j21
       do i = 1,i21
           p(i,j,1) = p(i,j,2)
           p(i,j,k20) = p(i,j,k19)
       enddo
       enddo
       !*** gauss-seidel method

       do k = 2,k19
       do j = 2,j19
       do i = 2,i19
           pcor = dijk*((p(i+1,j,k)+p(i-1,j,k))/(dx*dx) &
                       +(p(i,j+1,k)+p(i,j-1,k))/(dy*dy) &
                       +(p(i,j,k+1)+p(i,j,k-1))/(dz*dz)-q(i,j,k)) &
                       -p(i,j,k)
           g1 = g1 + pcor*pcor
           p(i,j,k) = p(i,j,k) + pcor
       enddo
       enddo
       enddo

       if( g1.le..001) then
           write(*,*)"kk=",kk
           exit
       endif
   enddo



[MPI実装例]

  do kk = 1,km
      g1 = 0.
      g_part=0.

      !*** boundary gondi'iioh for pressure
      !*** y-z
      do k = 1,k21
      do j = 1,j21
          p(1,j,k) = p(2,j,k)
          p(i20,j,k) = p(i19,j,k)
      enddo
      enddo
      !***x-z
      do k = 1,k21
      do i = 1,i21
          p(i,1,k) = p(i,2,k)
          p(i,j20,k) = p(i,j19,k)
      enddo
      enddo
      !*** x-y
      do j = 1,j21
      do i = 1,i21
          p(i,j,1) = p(i,j,2)
          p(i,j,k20) = p(i,j,k19)
      enddo
      enddo

      do ii=0,1

          if(ii .eq. 0) then

              !**********************!
              !! send to right side !!
              !**********************!
              if(myrank .ne. nproc-1) then
                  do itr_k=1, k21
                  do itr_j=1, j21
                      sendbuf2D(itr_j, itr_k) = p(iend(myrank), itr_j, itr_k)
                  enddo
                  enddo

                  call mpi_send(sendbuf2D, j21*k21, MPI_REAL, iup, 3, mpi_comm_world,ierr)
              endif

              if(myrank .ne. 0) then
                  call mpi_recv(recvbuf2D, j21*k21, MPI_REAL, idown, 3, mpi_comm_world,istatus,ierr)

                  do itr_k=1, k21
                  do itr_j=1, j21
                      p(istt(myrank)-1, itr_j, itr_k) = recvbuf2D(itr_j, itr_k)
                  enddo
                  enddo

              endif

          else

              !*********************!
              !! send to left side !!
              !*********************!
              if(myrank .ne. 0) then
                  do itr_k=1, k21
                  do itr_j=1, j21
                      sendbuf2D(itr_j, itr_k) = p(istt(myrank), itr_j, itr_k)
                  enddo
                  enddo
                  call mpi_send(sendbuf2D, j21*k21, MPI_REAL, idown, 4, mpi_comm_world, ierr)
              endif

              if(myrank .ne. nproc-1) then
                  call mpi_recv(recvbuf2D, j21*k21, MPI_REAL, iup, 4, mpi_comm_world,istatus,ierr)

                  do itr_k=1, k21
                  do itr_j=1, j21
                      p(iend(myrank)+1,itr_j,itr_k) = recvbuf2D(itr_j, itr_k)
                  enddo
                  enddo
              endif

          endif

          !****************************!
          !*** END OF DATA EXCHANGE ***!
          !****************************!



          !*** gauss-seidel method
          do k = 2,k19
          do j = 2,j19
          do i = istt(myrank)+ii ,iend(myrank) , 2   
              pcor = dijk*((p(i+1,j,k)+p(i-1,j,k))/(dx*dx) &
                          +(p(i,j+1,k)+p(i,j-1,k))/(dy*dy) &
                          +(p(i,j,k+1)+p(i,j,k-1))/(dz*dz)-q(i,j,k)) &
                          -p(i,j,k)
              g_part = g_part + pcor*pcor
              p(i,j,k) = p(i,j,k) + pcor
          enddo
          enddo
          enddo


      enddo
      call mpi_allreduce(g_part, g1, 1, MPI_REAL, MPI_SUM, MPI_COMM_WORLD, ierr)

      if( g1.le..001)then
          if(myrank.eq.0) write(*,'(a,I0)') ">> exit G-S loop coumnt=",kk
          exit
      endif

  enddo


変数iiにより、偶数奇数を判別しています。

反復計算の先頭にて、データの送受信が必要ですが、

  • 偶数ターンの時には、偶数番目から奇数番目へ、一つ上のiupへ、境界データを送信します。
  • 奇数ターンの時には、奇数番目から偶数番目へ、一つ下のidownへ、境界データを送信します。

データの送受信が済んだら、偶数ターンと奇数ターンでそれぞれ、ガウスザイデル法の計算を適用します。


[参考文献]

基礎からの流体力学! - 河村哲也著 3Dキャビティフロー

http://www.hpc.co.jp/news/images/seminar/3dcavity_problem.f90


この記事へのコメントをお寄せください

  • サイトへの書き込みに差し支えございましたら トータルディスクロージャーサイトサポート係へメールをお送りください
  • トータル・ディスクロージャ・サイトに投稿された文章と画像は、すべてその著作権がHPCシステムズ株式会社に帰属し、HPCシステムズ株式会社が著作権を所有することに同意してください。
  • あなたの文章が他人によって自由に編集、配布されることを望まない場合は、投稿を控えてください。
  • コメントを書き込む場合は名前にひらがなを織り交ぜてください。
  • あなたの投稿する文章と画像はあなた自身によって書かれたものであるか、パブリック・ドメインかそれに類する自由なリソースからの複製であることを約束してください。あなたが著作権を保持していない作品を許諾なしに投稿してはいけません!

<comments hideform="false" />


Comments

ノート:メニーコアと格闘ブログ:ポアソン方程式の並列化 その3

個人用ツール