OpenACCで2次元圧縮性流体コードをGPU対応させる

# 導入
この記事は2023年1月17日に国立天文台天文シミュレーションプロジェクト(CfCA)で行われるGPU講習会で使用するために書かれています。CfCAでは現在32機のNVIDIA A100を運用しており、天文シミュレーションなどに使用することが可能です。

以下、リンク先はCfCA登録者限定。
- [GPU講習会](https://www.cfca.nao.ac.jp/content/gpu_workshop_2022_2)
- [GPU使用申請](https://www.cfca.nao.ac.jp/content/gpu-%E5%8B%9F%E9%9B%86%E8%A6%81%E9%A0%85%E3%81%8A%E3%82%88%E3%81%B3%E5%88%A9%E7%94%A8%E7%94%B3%E8%AB%8B%E6%9B%B8)

この記事ではFortran90で書かれた2次元圧縮性流体コードをOpenACCを用いてGPU対応させます。サンプルコードは[こちら](https://github.com/cfcanaoj/gpuhydro/tree/master/KHf90openacc)です。以下では解説のためコードを簡単化していることがありますので、実際に動くコードをみたい場合は上記のサンプルを参考にしてください。またハンズオンで実際に作業したい方は[こちら](https://github.com/cfcanaoj/gpuhydro/blob/ws2022/KHf90openacc/PRACTICE.md)を見てください。

OpenACCを使うと以下のように`!$acc`から始まる指示行を用いて簡単にGPUに対応したコードを生成することができます。
~~~
!$acc 
~~~

# 方針
基本的に計算は全てGPU側（デバイス）で行い、ファイルの出力だけCPU側（ホスト）で行います。プログラムmainの中でサブルーチンは以下のように呼ばれます。グリッド生成`GenerateGrid`, でグリッドのデータをホストからデバイスに転送します。
また、問題生成, `GenerateProblem`, で流体変数のデータをホストからデバイスに転送します。`TimestepControl`で`dt`をデバイスからホストに `Output`で流体の変数をデバイスからホストに転送します。

~~~Fortran:
      call GenerateGrid
      call GenerateProblem
      call ConsvVariable
      do nhy=1,nhymax
         call TimestepControl
         call BoundaryCondition
         call StateVevtor
         call NumericalFlux1
         call NumericalFlux2
         call UpdateConsv
         call PrimVariable
         time=time+dt
         call Output
      enddo
~~~

# ホスト、デバイス間で共有するデータの宣言
このプログラムでは多くのデータをデバイス、ホスト間で共有します。例えば流体の変数`d,et,mv1,mv2,mv3`(それぞれ密度、全エネルギー密度、3方向の運動量密度)と断熱指数(`gam`)をデバイス、ホスト間で共有するには以下のように書きます。変数を宣言したあと`!$acc declare create (変数)`とします。
~~~Fortran:
      module modbasic
      implicit none
      real(8),dimension(in,jn,kn)::d,et,mv1,mv2,mv3

      real(8),parameter::gam=5.0d0/3.0d0
!$acc declare create(d,et,mv1,mv2,mv3)
!$acc declare create(gam)
      end module modbasic
~~~

このコードでは`in,jn,kn`は定数ですが、動的に割り付けるときもOpenACCの書き方は同じです。以下のようにしてください。
~~~Fortran:
      real(8),dimension(:,:,:),allocatable::d,et,mv1,mv2,mv3
!$acc declare create(d,et,mv1,mv2,mv3)
~~~
その後、ホストでallocateして`update device`(この命令が何かは後述)するとデバイス側で配列が正しく確保されます。
~~~Fortran:
    allocate(d(in,jn,kn))
!$acc update device (d)
~~~

# ホストからデバイスへの転送
こうしたデータをホストからデバイスに転送するのは簡単です。以下のようにします。以下のコードではdo loopはホストで実行され、密度`d`　の値が更新されます。`!$acc update device (d)`で密度`d`が転送されます。

~~~Fortran:
      subroutine GenerateProblem
      use modbasic
      implicit none
      do k=ks,ke
      do j=js,je
      do i=is,ie
          d(i,j,k) = 1.0d0
      enddo
      enddo
      enddo
!$acc update device (d)
      return
      end subroutine GenerateProblem
~~~

# デバイスからホストへの転送
逆にデバイスからホストに転送するのも簡単です。以下のようにします。`Output`では`time`が 前回ファイルをアウトプットした時間`tout`から`dtout`分経過したときにファイルを出力します。

その場合、`return`の後の行の`!$acc update host (d,v1,v2,v3,p)`が呼ばれます。この時密度や運動量などの変数の値がホスト側に転送されるので。その値をファイルに出力しています。

~~~Fortran:
      subroutine Output
      use modbasic
      implicit none
      character(40)::filename
      real(8),save::tout
      data tout / 0.0d0 / 
      real(8),parameter:: dtout=1.0d-2
      integer::nout
      data nout / 1 /
      if(time .lt. tout+dtout) return
!$acc update host (d,v1,v2,v3,p)

      write(filename,'(a2,i5.5,a4)')"Sc",nout,".xss"
      filename = trim(dirname)//filename
 open(unitout,file=filename,status='replace',form='formatted') 
      write(unitout,*) "# ",time
      k=ks
      do j=js,je
      do i=is,ie
         write(unitout,'(7(1x,E12.3))') x1b(i),x2b(j),d(i,j,k),v1(i,j,k),v2(i,j,k),v3(i,j,k),p(i,j,k)
      enddo
         write(unitout,*)
      enddo
      close(unitout)

      write(6,*) "output:",nout,time

      nout=nout+1
      tout=time

      end subroutine Output
~~~
## Unified Memory
このデータ輸送部分は結構面倒です。ひとまずloopの変更に集中したい場合は`-gpu=managed`というコンパイルオプションをつけるとデータの輸送が自動的に行われます。ただし、処理速度が出ないことが多いので最終的には外すことになるでしょう。

# ループの独立実行

## ループの独立実行
いよいよGPUによる計算パートの書き方を説明します。今、`d, nflux1, nflux2, x1a, x2a, dt`などは全てデバイス側で持っているとします。

まず、デバイス側で実行したいパートの指定には`!$acc kernels`と`!$acc end kernels`を用います。これで囲んだパートはデバイス側で実行されます。

下記の例ではdo loopがあります。`!$acc loop independent`をループの前に書くと、そのループについては並列に実行しても良いという指示になります。この場合3ループありますが、`collapse(3)`とするとこで3ループを結合してくれます。collapseを明示しなくてもコンパイラが自動でループをまとめてくれることもありますが、明示的につけたほうが安全です。もし`collapse`を使わない場合はループ一つ一つに`!$acc loop independent`をつけてください。

注意するべきなのは、このループの中に出てくる変数はスカラーと`i,j,k`がついている配列だけということです。この場合は`!$acc loop independent`だけで並列実行してくれます。

~~~Fortran:
      subroutine UpdateConsv
      use modbasic
      use fluxmod
      implicit none
      integer::i,j,k

!$acc kernels
!$acc loop collapse(3) independent
      do k=ks,ke
      do j=js,je
      do i=is,ie
         
         d(i,j,k) = d(i,j,k)  &
     & +dt*( &
     &  (- nflux1(mden,i+1,j,k) &
     &   + nflux1(mden,i  ,j,k))/(x1a(i+1)-x1a(i))  &
     & +(- nflux2(mden,i,j+1,k) &
     &   + nflux2(mden,i,j  ,k))/(x2a(j+1)-x2a(j))  &
     &      )

      enddo
      enddo
      enddo
!$acc end kernels

      return
      end subroutine UpdateConsv
~~~

GPUではループが並列的に実行されるので同期をとりたい区切りで`!$acc kernels`と`!$acc end kernels`を用いてください。つまり以下の例ではpartAが終わった時に全ての`A`がただしく更新されていますが、
~~~Fortran:
!============
! part A
!============
!$acc kernels
!$acc loop independent
      do i=1,n
        A(i) = A(i)+ B(i)
      enddo
!$acc end kernels
!============
! part B
!============
!$acc loop independent
      do i=1,n
        C(i) = A(i)+ B(i)
      enddo
!$acc end kernels
~~~
以下の例では最初のループを計算し終わったスレッドから次々に二つ目のループを計算しにいくので`A`が正しくアップデートされていない可能性があります。
~~~Fortran:
!$acc kernels
!$acc loop independent
      do i=1,n
        A(i) = A(i)+ B(i)
      enddo
!$acc loop independent
      do i=1,n
        C(i) = A(i)+ B(i)
      enddo
!$acc end kernels
~~~

## ループの独立実行（private配列入り）
やや難しいケースとして`i,j,k`がついてない配列がある場合について説明します。以下のサブルーチンはセルセンターにある変数を用いてセルエッジでの変数の値を内挿によって求めています。

ここで気をつけて欲しいのは配列`dsv(:), dsvp(:), dsvm(:)`には`i,j,k`がついていませんが、意味的には明らかに`i,j,k`について独立にそれぞれ計算されるものであるということです。この変数は並列実行のそれぞれのthread間で共有されていないprivateな変数であると教えてあげる必要があります。そこで`private(dsv,dsvp,dsvm)`という指示がされています。

~~~Fortran:
      subroutine NumericalFlux1
      use modbasic, only: is,ie,in,js,je,jn,ks,ke,kn,gam
      use fluxmod
      implicit none
      integer::i,j,k
      real(8),dimension(nhyd):: dsvp,dsvm,dsvc,dsv
      real(8),dimension(nhyd,in,jn,kn):: leftpr,rigtpr
      real(8),dimension(2*mflx+madd,in,jn,kn):: leftco,rigtco
      real(8),dimension(2*mflx+madd):: leftst,rigtst
!$acc kernels
      k=ks
!$acc loop collapse(2) independent private(dsv,dsvp,dsvm)
      do j=js,je
      do i=is-1,ie+1
         dsvp(:) = (svc(:,i+1,j,k) -svc(:,i,j,k)                 )
         dsvm(:) = (                svc(:,i,j,k) - svc(:,i-1,j,k))

         call vanLeer(dsvp,dsvm,dsv)
         leftpr(:,i+1,j,k) = svc(:,i,j,k) + 0.5d0*dsv(:)
         rigtpr(:,i  ,j,k) = svc(:,i,j,k) - 0.5d0*dsv(:)
      enddo
      enddo
!$acc end kernels
      return
      end subroutine Numericalflux1
~~~

collapseを明示しなくてもコンパイラが自動でループをまとめてくれることもありますが、private配列がある場合はまとめてはくれません。つまり以下のようにまとめて書いてはいけません。
~~~Fortran:
!$acc loop independent private(dsv,dsvp,dsvm)
      do j=js,je
      do i=is-1,ie+1
~~~

このプログラムの中ではもう1箇所（2次元なので実際は2倍）プライベート変数が指定されている箇所があります。セルエッジの右側の物理量と左側の物理量がわかっているときに数値フラックスを求める部分です。こういった箇所では少し気を付ける必要があります。
~~~Fortran:
      subroutine NumericalFlux1
      use modbasic, only: is,ie,in,js,je,jn,ks,ke,kn,gam
      use fluxmod
      implicit none
      integer::i,j,k
!$acc kernels
!$acc loop collapse(2) independent private(leftst,rigtst,nflux)
      do j=js,je
      do i=is,ie+1
         leftst(:)=leftco(:,i,j,k)
         rigtst(:)=rigtco(:,i,j,k)
         call HLLE(leftst,rigtst,nflux)
         nflux1(mden,i,j,k)=nflux(mden)
         nflux1(mrv1,i,j,k)=nflux(mrvu)
         nflux1(mrv2,i,j,k)=nflux(mrvv)
         nflux1(mrv3,i,j,k)=nflux(mrvw)
         nflux1(meto,i,j,k)=nflux(meto)
      enddo
      enddo
!$acc end kernels

      return
      end subroutine Numericalflux1
~~~

## ループの独立実行の中での関数コール
先ほどの数値フラックスを求める例では並列実行されているループの中で関数コール`HLLE`がされております。この部分をデバイス側で実行するためには以下のように`module`の`contains`の中で関数を定義して関数名のあとに`!$acc routine seq`と書く必要があります。これはデバイス側で並列実行ではなくて逐次実行(sequential)されるという指示行です（この関数の中でさらに並列にすることもできますが、ここでは触れません）。呼び出し側が同じファイル内の場合はmoduleで囲まなくても良いようです。

~~~Fortran:
      module numfluxmod
      contains 
      subroutine HLLE(leftst,rigtst,nflux)
!$acc routine seq
      use fluxmod
      implicit none
      real(8),dimension(2*mflx+madd),intent(in)::leftst,rigtst
      real(8),dimension(mflx),intent(out)::nflux
      real(8),dimension(mflx)::ul,ur,fl,fr
      real(8)::csl,csr
      real(8):: vl, vr
      real(8):: sl, sr

      ul(1:mflx) = leftst(1:mflx)
      fl(1:mflx) = leftst(mflx+1:2*mflx)
      ur(1:mflx) = rigtst(1:mflx)
      fr(1:mflx) = rigtst(mflx+1:2*mflx)
      csl=leftst(mcsp)
      csr=rigtst(mcsp)
       vl=leftst(mvel)
       vr=rigtst(mvel)

       sl = min(vl,vr) - max(csl,csr)
       sl = min(0.0d0,sl)
       sr = max(vl,vr) + max(csl,csr)
       sr = max(0.0d0,sr)

       nflux(:) = (sr*fl(:)-sl*fr(:) +sl*sr*(ur(:)-ul(:)))/(sr-sl)

      return
      end subroutine HLLE
      end module numfluxmod
~~~

このHLLEは元々`!$acc kernels`-`!$acc end kernels`内で呼ばれていることに気をつけてください。つまり以下のような構造になっています。
~~~Fortran:    
    subroutine RoutineA
!$acc kernels
!$acc loop independent 
      do 
         call RoutineB
      enddo
!$acc end kernels
    end RoutineA

      subroutine RoutineB
!$acc routine seq
      end subroutine RoutineB
~~~

# ループのリダクション
前節のようにループを独立に実行するのは比較的簡単ですが、合計、最大、最小を求めるなどのリダクションを並列実行するのは一般的には難しいです。とはいえ、OpenACCでのリダクションはCudaに比べれば圧倒的に簡単です。

このコードの場合、次のタイムステップを決めるために音速＋流速が最も速いところを見つける必要があります。`!$acc loop collapse(3) reduction(min:dtmin) `をループの前につけることでそのリダクションを実行してくれます。ちなみに最大の場合は`reduction(max:変数)`で合計は`reduction(+:変数)`です。 ~~変数は配列は許されないので、配列を使っている場合は少し書き方を工夫する必要があります。~~ [OpenAcc 2.7からは配列も許されます](https://www.openacc.org/blog/whats-new-openacc-27)。

~~~Fortran:
      subroutine TimestepControl
      use modbasic
      implicit none
      real(8)::dtl1, dtl2, dtl3, dtlocal
      real(8)::dtmin
      integer::i,j,k
!$acc kernels
      dtmin=1.0d90
!$acc loop collapse(3) reduction(min:dtmin)      
      do k=ks,ke
      do j=js,je
      do i=is,ie
         dtl1 =(x1a(i+1)-x1a(i))/(abs(v1(i,j,k)) +cs(i,j,k))
         dtl2 =(x2a(j+1)-x2a(j))/(abs(v2(i,j,k)) +cs(i,j,k))
         dtlocal = min (dtl1,dtl2)
         if(dtlocal .lt. dtmin) dtmin = dtlocal
      enddo
      enddo
      enddo   
      dt = 0.05d0 * dtmin
!$acc end kernels
!$acc update host (dt)
      
      return
      end subroutine TimestepControl
~~~

僕がreductionで困ったことも紹介します。reductionにはモジュール変数を使わないほうが良いです。

~~~Fortran:
module varmod
  implicit none
  real(8):: Amod
!$acc declare create(Amod)
end module varmod

program main
  use varmod
  implicit none
  integer,parameter:: nmax=3
  integer,parameter:: imax=100
  integer:: i,n
  real(8):: Aloc
!$acc declare create(Aloc)

  
  write(6,*) "TEST1"
  do n=1,nmax
!$acc kernels
     Aloc = 0.0d0
!$acc loop reduction(+:Aloc)  
     do i=1,imax
        Aloc = Aloc + 1.0d0
     enddo
!$acc end kernels
!$acc update host (Aloc)
     write(6,*) "iter=",n," sum=",Aloc
  enddo

  write(6,*) "TEST2"
  do n=1,nmax
!$acc kernels
     Amod = 0.0d0
!$acc loop reduction(+:Amod)
     do i=1,imax
        Amod = Amod + 1.0d0
     enddo
!$acc end kernels
!$acc update host (Amod)
     write(6,*) "iter=",n," sum=",Amod
  enddo

  write(6,*) "TEST3"
  do n=1,nmax
     Amod = 0.0d0
!$acc update device (Amod)
!$acc kernels
!$acc loop reduction(+:Amod)
     do i=1,imax
        Amod = Amod + 1.0d0
     enddo
!$acc end kernels
!$acc update host (Amod)
     write(6,*) "iter=",n," sum=",Amod
  enddo
  stop
end program main
~~~
上記のコードを実行すると実は以下のようになります。モジュール変数`Amod`を使ったTEST2の結果が思ったものではないですね。モジュール変数をリダクションで使うと、`Amod=0`が反映されず、前のイテレーションの`Amod`が初期値になってしまいます。なのでreducitonにはローカル変数を使うほうがよいですし、もしモジュール変数を使う場合にはTEST3のようにホストで初期化してそれをデバイスでアップデートしてください。
~~~
 TEST1
 iter=            1  sum=    100.0000000000000     
 iter=            2  sum=    100.0000000000000     
 iter=            3  sum=    100.0000000000000 

 TEST2
 iter=            1  sum=    100.0000000000000     
 iter=            2  sum=    200.0000000000000     
 iter=            3  sum=    300.0000000000000

 TEST3
 iter=            1  sum=    100.0000000000000     
 iter=            2  sum=    100.0000000000000     
 iter=            3  sum=    100.0000000000000
~~~
# コンパイルと実行
CfCAの環境ではコンパイラは以下のように指定します。
~~~make:Makefile
exe=kh.x
fc= nvfortran
foptopenacc = -Minfo=accel -acc
fopt = -g -traceback -O2

${exe}: main.f90
	${fc} ${fopt} ${foptopenacc} $< -o ${exe}
~~~
問題サイズが大きくなって（2GB以上の固定配列を持とうとすると？）エラーが出るときには、以下の`-mcmodel=medium`を指定する。
~~~make:Makefile    
fc= nvfortran -mcmodel=medium
~~~

CfCAのGPUサーバ群ではslurmを利用していますので、コンパイルが出来たら以下のように実行します。バッチファイルの中身は[こちら](https://github.com/cfcanaoj/gpuhydro/blob/master/KHf90openacc/sj_g00.sh)を見てください。
~~~bash
sbatch sj_g00.sh
~~~

# 計測
GPUコードを書いたからには計測してみたくなります。A100環境とSkylake 2CPUを比較してみました。計測にはopenmpの関数を用いました。

~~~Fortran:
      program main
      use omp_lib
      implicit none
      real(8)::time_begin,time_end
! main loop
      time_begin = omp_get_wtime()
      do nhy=1,nhymax
         call TimestepControl
         call BoundaryCondition
         call StateVevtor
         call NumericalFlux1
         call NumericalFlux2
         call UpdateConsv
         call PrimVariable
         time=time+dt
      enddo
      time_end = omp_get_wtime()
      write(6,*) "sim time [s]:", time_end-time_begin
      write(6,*) "time/count/cell", (time_end-time_begin)/(ngrid**2)/nhymax
      end program main
~~~

下記は深く考えず計っただけなので、あとで修正するかもしれませんが、以下のようになりました。Skylakeではopenmpを用いて並列化しています（[コードはこちら](https://github.com/cfcanaoj/gpuhydro/tree/master/KHf90openmp)）。

|問題サイズ|2 Skylake | A100|
| --- | --- | --- |
|500x500| 170 s | 21.9 s|

メモリバンド幅はSkylakex2が256 GB/sでA100が2000 GB/sなのでその比が直接表れているかもしれないです。問題サイズは500^2で20000stepの計算です。


# MPIを利用する。
さらに大規模にジョブを走らせるためにはMPI化する必要があります。その場合、以下の二点に気をつけてコードを変更する必要があります。サンプルコードは[こちら](https://github.com/cfcanaoj/gpuhydro/tree/mpi/KHf90openacc)です。上記のものの違うブランチのものとなっています。

## 1 MPI processに1GPUを対応させる
下記のように`acc_set_device_num`を使ってプロセスにGPUのデバイスを割り当てます。
~~~Fortran:
  use openacc

  ngpus = acc_get_num_devices(acc_device_nvidia)
  if(myid_w == 0) then
     print *, "num of GPUs = ", ngpus
  end if

  gpuid = mod(myid_w, ngpus)
  if(ngpus == 0) gpuid = -1
  if(gpuid >= 0) then
     call acc_set_device_num(gpuid, acc_device_nvidia)
  end if
~~~

## MPI関数を使うときにdeviceのメモリを指定する。
上記のモデルではデバイス側に全てのデータがあると想定しています。デバイス側のメモリにある変数をMPI変数で受け渡すには以下のようにします。

~~~Fortran:
!$acc host_data use_device(varsendXstt,varrecvXst)
     nreq = nreq + 1         
     call MPI_IRECV(varrecvXstt,mgn*jn*kn*nbc &
    & , MPI_DOUBLE &
    & , n1m,1100, comm3d, req(nreq), ierr)

     nreq = nreq + 1
     call MPI_ISEND(varsendXstt,mgn*jn*kn*nbc &
    & , MPI_DOUBLE &
    & , n1m, 1200, comm3d, req(nreq), ierr)

     if(nreq .ne. 0) call MPI_WAITALL ( nreq, req, stat, ierr )
     nreq = 0
!$acc end host_data
~~~
## バッチジョブで使用するGPUの枚数を指定する
CfCAの場合、以下のようにバッチスクリプトを書きます。以下では4MPI processでGPU4枚使用する場合です。プロセス数を変えたいときは4を書き換えてください。
~~~Slurm:
#SBATCH --ntasks=4
#SBATCH --gres=gpu:4

module load nvhpc

mpiexec -n 4 ./kh.x > log.dat
~~~

## MPI版の計測
GPU4枚を使った場合とGPU1枚を使った場合の計算時間の比較を以下の表にまとめました。スピードアップは計算時間の比で理想的には4倍になります。しかし、500^2の問題サイズをみると1.9倍にしかスピードアップされていません。2次元の問題だけあって、問題サイズが小さめであり、通信速度が無視できないようです。4GPUの計算力を生かし切るには2000^2以上の問題サイズが必要なようです。

|問題サイズ|A100|A100x4|スピードアップ|
| --- | --- | --- |--- |
| 500^2|21 s  | 11 s  | 1.9  |
|1000^2|76.2 s|	26.8 s|	2.84|
|2000^2	|306 s|	83.2 s|	3.67|



# 参考リンク
- [第1回 High Performance Computing Physics (HPC-Phys) 勉強会](https://hpc-phys.kek.jp/workshop/workshop180822.html)
- [第2回 High Performance Computing Physics (HPC-Phys) 勉強会](https://hpc-phys.kek.jp/workshop/workshop181201.html)
- [東京大学第157回並列プログラミング講習会「GPUプログラミング入門」](https://www.cc.u-tokyo.ac.jp/events/lectures/157/)
- [東京大学第167回並列プログラミング講習会「OpenACCとMPIによる
マルチGPUプログラミング入門」](https://www.cc.u-tokyo.ac.jp/events/lectures/167/)