導入

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

以下、リンク先はCfCA登録者限定。

この記事ではFortran90で書かれた2次元圧縮性流体コードをOpenACCを用いてGPU対応させます。サンプルコードはこちらです。以下では解説のためコードを簡単化していることがありますので、実際に動くコードをみたい場合は上記のサンプルを参考にしてください。またハンズオンで実際に作業したい方はこちらを見てください。

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

!$acc

方針

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

      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

ホスト、デバイス間で共有するデータの宣言

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

      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の書き方は同じです。以下のようにしてください。

      real(8),dimension(:,:,:),allocatable::d,et,mv1,mv2,mv3
!$acc declare create(d,et,mv1,mv2,mv3)

その後、ホストでallocateしてupdate device(この命令が何かは後述)するとデバイス側で配列が正しく確保されます。

    allocate(d(in,jn,kn))
!$acc update device (d)

ホストからデバイスへの転送

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

      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)が呼ばれます。この時密度や運動量などの変数の値がホスト側に転送されるので。その値をファイルに出力しています。

      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だけで並列実行してくれます。

      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がただしく更新されていますが、

!============
! 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が正しくアップデートされていない可能性があります。

!$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)という指示がされています。

      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配列がある場合はまとめてはくれません。つまり以下のようにまとめて書いてはいけません。

!$acc loop independent private(dsv,dsvp,dsvm)
      do j=js,je
      do i=is-1,ie+1

このプログラムの中ではもう1箇所（2次元なので実際は2倍）プライベート変数が指定されている箇所があります。セルエッジの右側の物理量と左側の物理量がわかっているときに数値フラックスを求める部分です。こういった箇所では少し気を付ける必要があります。

      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で囲まなくても良いようです。

      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内で呼ばれていることに気をつけてください。つまり以下のような構造になっています。

    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からは配列も許されます。

      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にはモジュール変数を使わないほうが良いです。

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の環境ではコンパイラは以下のように指定します。

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を指定する。

Makefile

fc= nvfortran -mcmodel=medium

CfCAのGPUサーバ群ではslurmを利用していますので、コンパイルが出来たら以下のように実行します。バッチファイルの中身はこちらを見てください。

sbatch sj_g00.sh

計測

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

      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を用いて並列化しています（コードはこちら）。

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

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

MPIを利用する。

さらに大規模にジョブを走らせるためにはMPI化する必要があります。その場合、以下の二点に気をつけてコードを変更する必要があります。サンプルコードはこちらです。上記のものの違うブランチのものとなっています。

1 MPI processに1GPUを対応させる

下記のようにacc_set_device_numを使ってプロセスにGPUのデバイスを割り当てます。

  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変数で受け渡すには以下のようにします。

!$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を書き換えてください。

#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

OpenMPを使いたい

OpenACCのディレクティブをOpenMPのディレクティブに書き換えることができます。

ツールを利用することもでき、便利です。2025年7月30日現在、下記のローカル変数のdeclareのところや

deviceで処理する場所の指定。

!$acc kernel
   .....
!$acc end kernel

!$omp target
   .....
!$omp end target

変数宣言(グローバル変数、module内)：

    module testmod
    real(8) a
!$acc declare create (a)
    end module testmod

    module testmod
    real(8) a
!$omp declare target (a)
    end module testmod

変数宣言(ローカル変数、subroutne内)：

    subroutine testsub
    real(8) a
!$acc data create (a)
    ...
!$acc end data
    end subroutine testsub

    subroutine testsub
    real(8) a
!$omp target data map(alloc:a)
    ...
!$omp end target data
    end subroutine testsub

変数更新：

!$acc update device (a)

!$omp update target to (a)

ループ：
３次元計算で使いそうな3重ループの場合。

!$acc loop collapse(3) independent

!$omp loop order(concurrent) collapse(3)

!$acc loop collapse(3) reduction(+:var)

!$omp target teams distribute parallel do collapse(3) reduction(+:var)

逐次実行する関数の呼び出し：

!$acc routine seq
subroutine myfunc(x)
...
end subroutine

!$omp declare target 
subroutine myfunc(x)
...
end subroutine
!$omp end declare target

MPI等の外部ルーチンへ変数を受け渡す：
３次元計算で使いそうな3重ループの場合。

!$acc host_data use_device(a)
call MPI_IRECV(a,...)
!$acc end host_data

!$omp target data map (to:a) use_device_ptr(a)
call MPI_IRECV(a,...)
!$omp end target data

参考リンク