導入
この記事は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
ホスト、デバイス間で共有するデータの宣言
このプログラムでは多くのデータをデバイス、ホスト間で共有します。例えば流体の変数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の環境ではコンパイラは以下のように指定します。
exe=kh.x
fc= nvfortran
foptopenacc = -Minfo=accel -acc
fopt = -g -traceback -O2
${exe}: main.f90
${fc} ${fopt} ${foptopenacc} $< -o ${exe}
問題サイズが大きくなって（2GB以上の固定配列を持とうとすると？）エラーが出るときには、以下の-mcmodel=mediumを指定する。
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 |