導入

いままでいくつか圧縮性流体の方程式をメッシュ法で解くコードのGPU化について記事を書いてきました。

GPU対応の難しいところは、使うハードウェアの製造元によってプログラミング言語を変えなくてはいけないことです。汎用で使える言語もありますが、使用例が少なく一般に習得難易度が高くなってしまっています。

僕は普段、Fortranでコードを書いていまして、その場合一番簡単なのはNVIDIAのGPUを使い、Fortran + OpenACCを使うことです。一つ目の記事はそれを意識したものとなっています。

一方でAMDのGPUを使いたい場合もあるかと思います。その時にはOpenACCが使えないので代わりにOpenMPを使用してみたけれども、コンパイラが未成熟で今のところ難しかったというのが二つ目の記事です。ただしこの記事の実行環境はNVIDIAのGPU, コンパイラです。AMDのものだとどうなるのかは未確認です。

OpenMPを使用した場合でも、Ｃ言語の場合はFortranよりはよく最適化してくれるみたいなので、三つ目の記事ではC++を使ってコードを書いています。

本稿では、三つ目の記事を発展させて、C++とOpenMPの枠組みでMPIも使って2次精度の磁気流体コードを書きます。数値フラックスにはHLLDを使用しており、div B=0を成立させるためにdivergence cleaning法を用いています。

詳細は後で紹介しますが、大きく分けて二つ困難がありました。一つは計算をGPU device上で完結させたいのにOpenMPのコードはhostと通信したがるということです。もう一つはC++の特徴であるクラスや構造体を導入するとGPU対応させるときの書き方が難しくなるというものでした。

計算コードはこちらで公開しています。現状以下のよう10のコードを公開していますが、MPI化もされているもの、３Ｄのものが最新です。参考にするならそちらをご覧ください。他のものは簡単のためや理解のため、あえてMPI化されてないものが見たい場合にご使用ください。

Problem	Language	parallel	Name
2D Kelvin-Helmholtz Instability	Fortran	OpenMP(CPU)	KHf90openmp
2D Kelvin-Helmholtz Instability	Fortran	MPI, OpenMP(CPU)	KHf90openmp_mpi
2D Kelvin-Helmholtz Instability	Fortran	OpenACC	KHf90openacc
2D Kelvin-Helmholtz Instability	Fortran	MPI, OpenACC	KHf90openacc_mpi
3D upwind advection	C++	OpenMP(GPU)	ADcppopenmp
3D upwind advection	C++	OpenACC	ADcppopenacc
3D magneto-hydrodynamic deacaying turbulence	Fortran	OpenACC	DTf90openacc
3D magneto-hydrodynamic deacaying turbulence	Fortran	MPI, OpenACC	DTf90openacc_mpi
3D magneto-hydrodynamic deacaying turbulence	C++	OpenMP(GPU)	DTcppopenmp
3D magneto-hydrodynamic deacaying turbulence	C++	MPI, OpenMP(GPU)	DTcppopenmp_mpi

実装の解説

以下ではDTcppopenmp_mpiの実装の解説をしていきます。C++で書かれた３次元の磁気流体コードでMPI, GPU向けOpenMPで並列化されています。

コンパイル　オプション

NVIDIAのコンパイラを使っていまして、コンパイルオプションは以下のようになっています。-mp=gpuでOpenMPを使ってGPU対応をします。オプションの-Minfo=mp,accelが大事です。これをつけているとGPU offloadしたループごとにhost, device間でどのように通信しているなどの情報がでます。-gpu=mem:separateはhost-device間の通信を手動でするというものです。とはいえ後述するように、コンパイラが勝手に通信を挟んできます。

coptgpu := -mp=gpu -Minfo=mp,accel -gpu=cc80,mem:separate

GPU対応でのデータの取り扱いの一般論

OpenMPはデータを(1)enter allocして(2)update toして(3)update fromして(4)exit deleteするのが基本セットと覚えておきましょう。メモリの解放は省略しても良いです。namespace内の変数の場合は(0)#pragma omp declare targetと#pragma omp end declare targetで囲みましょう。

// (0) namespace内の変数の場合のみ　そうでなければ不要
#pragma omp declare target
    double *a;
#pragma omp end declare target

// (1) メモリ確保
    a = (double*) malloc(sizeof(double)*num);
#pragma omp target enter data map(alloc: a[0:num])

// (2) CPU → GPU のデータ転送
    for (int n=0;n<num;n++){a[n] = 1.0e0;}
#pragma omp target update to(a[0:num])

// (3) GPU → CPU のデータ転送
    for (int n=0;n<num;n++){a[n] = 2.0e0;}
#pragma omp target update from(a[0:num])

// (4) メモリ解放
#pragma omp target exit data map(delete:a[0:num])

多次元配列の取り扱い

C++では多次元配列をどう扱うのかというのが、ひとつ問題となります。実のところGPU化を念頭に置くと素朴に配列を1次元で確保してそのまま使うというのが一番簡単です。例えば3次元のin x jn x kn の配列にアクセスするとき以下のようにするということです。

    a[((k)*jn+j)*in+i] = 3.0e0; //(1)素朴な書き方

ただ、上記のものは見た目が悪く、よくある実装では以下の(3)のように構造体使って疑似的に多次元配列化のようにアクセスします。

    a.data[((k)*jn+j)*in+i] = 3.0e0;// (2) aというのクラスのメンバ変数dataにアクセス
    a(k,j,i)= 3.0e0; // (3) このように簡単に書ける

一方でGPU上のメモリを確保するときに結局生の配列の情報を書かないといけないので、構造体で隠蔽しているデータの構造を把握しておく必要がでてきて、(2)のようなものと理解しておく必要があるし、MPI関係でメンバ変数を使えない場面もでてくるので、書き方を工夫する必要がでてきます。その意味では(1)のように書いてしまうのも良いかもしれませんし、プリプロセッサを使って(k,j,i)を[((k)*jn+j)*in+i]に変換したりするのも視野に入りますが、苦労の末構造体を使ってもできたのでここではその結果を紹介することにします。

流体場の変数はprimitive variables用にP,conservative variables用に U, numerical flux用にFx,Fy,Fzを用意しています。

hydro.hpp

  extern FieldArray<double> P; //! P(nprim ,ktot,jtot,itot)
  extern FieldArray<double> U; //! U(mconsv,ktot,jtot,itot)
  extern FieldArray<double> Fx,Fy,Fz;

FieldArrayの定義は以下のようになっています。メンバー変数としてdata,nv,n3,n2,n1,sizeが存在しています。ここで各方向のグリッドサイズnv,n3,n2,n1と全サイズsizeはGPU用の指示行で陽にでてくるのでpublicとして定義しておくのが重要です。operatorの定義をしているので上記の(3)のようなアクセスが可能になります。

hydro.hpp

  template <typename T>
  class FieldArray {
  public:
    T* data = nullptr;
    int nv = 0, n3 = 0, n2 = 0, n1 = 0;
    size_t size = 0;
    
    void allocate(int _nv, int _n3, int _n2, int _n1) {
      nv = _nv; n3 = _n3; n2 = _n2; n1 = _n1;
      size = static_cast<size_t>(nv) * n3 * n2 * n1;
      data = static_cast<T*>(malloc(sizeof(T) * size));
  }
        
    inline const T& operator()(int n, int k, int j, int i) const noexcept {
      return data[((n*n3 + k)*n2 + j)*n1 + i];
    }
    inline  T& operator()(int n, int k, int j, int i)  noexcept {
      return data[((n*n3 + k)*n2 + j)*n1 + i];
    }
  };

この変数を以下のように確保します。Uに特化して説明します。Uはnamespace hydflux_modの中で定義されるglobal変数なので#pragma omp declare target-#pragma omp end declare targetで囲みます。

AllocateHydroVariablesの引数にUを指定しなくてもnamespace上で定義されているので、問題ないのですが、引数で指定しているとコンパイラが全体的な構造を理解しやすいようで、無駄なhost-device間の通信を避けられることがあるようです。CPU上でメモリを確保した後、GPU上でメモリを確保するために具体的に#pragma omp target enter data map (alloc: U.data[0: U.size])としているところに特徴があります。U.data[0: U.size]は一次元の配列であったことを思い出してください。コロンの前は最初のindex,コロンの後は全体のサイズです。初期化したあと#pragma omp target update to ( U.data[0: U.size], U.n1, U.n2, U.n3, U.nv)で配列のデータ数に関するメンバ変数も全てdevice側に伝えます。

hydro.cpp

namespace hydflux_mod {
#pragma omp declare target
  FieldArray<double> U; //! U(mconsv,ktot,jtot,itot)
#pragma omp end declare target

  void AllocateHydroVariables(FieldArray<double>& U){
  
  U.allocate(mconsv,ktot,jtot,itot);// CPU上でallocate
#pragma omp target enter data map (alloc: U.data[0: U.size]) //GPU上でallocate
  for (int m=0; m<mconsv; m++)
  for (int k=0; k<ktot; k++)
  for (int j=0; j<jtot; j++)
  for (int i=0; i<itot; i++) {
  	   U(m,k,j,i) = 0.0;
  }
#pragma omp target update to ( U.data[0: U.size], U.n1, U.n2, U.n3, U.nv)
  }

}

GPU offload

numerical fluxを計算するところではGPU offloadで高速化が可能です。GetNumericalFlux1の引数に流体の場の変数を書く必要はありませんが、書いたコードのほうが最適化されやすかったのでわざわざ書いています。
この関数内は3つのパートからできています。A-partでcell centerからcell edgeへprimitive variablesを補間しています。van Leer limiterもここでかけています。その後B-partでconservative variablesのU, Fを計算します。さらにC-partなのでHLLDを使ってnumerical fluxを計算します。

実はGPU化する前のプログラムでは上記のA-partの結果とB-partの結果を一度３次元のindexがついた配列に保存していたのですが、そういうglobalっぽい変数があるとhost-device間で無駄な通信が発生したこともあり、アルゴリズムを変更しました。なるべく関数内では新しくglobalっぽい変数を出さないほうが良いみたいです。

hydro.cpp

void GetNumericalFlux1(const GridArray<double>&G,const FieldArray<double>& P,FieldArray<double>& Fx){

#pragma omp target teams distribute parallel for collapse(2)
    for (int k=ks; k<=ke; k++)
    for (int j=js; j<=je; j++){
	  for (int i=is; i<=ie+1; i++) {
      // A-part: -----------------------------------
	  // Pick up related cell
	  double Pleftc1[nprim];
	  double Pleftc2[nprim];
	  double Prigtc1[nprim];
	  double Prigtc2[nprim];  
	  for (int n=0; n<nprim; n++){
	    Pleftc1[n] = P(n,k,j,i-2);
	    Pleftc2[n] = P(n,k,j,i-1);
	    Prigtc1[n] = P(n,k,j,i  );
	    Prigtc2[n] = P(n,k,j,i+1);
	  }
	  // Calculte Left state
	  double Plefte [nprim];
	  double dsvp,dsvm,dsv;
	  /* | Pleftc1   | Pleftc2 =>| Prigtc1   | Prigtc2   |  */
	  /*                     You are here                   */
	  for (int n=0; n<nprim; n++){
	    dsvp =  Prigtc1[n]- Pleftc2[n];
	    dsvm =              Pleftc2[n]- Pleftc1[n];
	    vanLeer(dsvp,dsvm,dsv);
	    Plefte[n] = Pleftc2[n] + 0.5e0*dsv;
	  }
      // B-part: -----------------------------------
	  double Clefte [2*mconsv+madd];
	Clefte[mudn] = Plefte[nden]; // rho
    ...
	Clefte[mfdn] = Plefte[nden]*Plefte[nve1];
    ...
    double Prigte [nprim];
	  for (int n=0; n<nprim; n++){
	    dsvp =  Prigtc2[n]- Prigtc1[n];
	    dsvm =              Prigtc1[n]- Pleftc2[n];
	    vanLeer(dsvp,dsvm,dsv);
	    Prigte[n] = Prigtc1[n] - 0.5e0*dsv;
	  }
	  double Crigte [2*mconsv+madd];

	Crigte[mudn] = Prigte[nden]; // rho
    ...
	Crigte[mfdn] = Prigte[nden]*Prigte[nve1];
      // C-part: -----------------------------------
	 double numflux [mconsv];
	 HLLD(Clefte,Crigte,numflux);
	 Fx(mden,k,j,i) = numflux[mden];
    ...
	}// i-loop
  }// j,k-loop
}

Reduction

陽解法の流体計算では時間ステップdtを計算するためにreductionが必要になります。reductionそのものは簡単で#pragma omp target teams distribute parallel for reduction(min:dtminl) collapse(3)で終わりです。ただし、ここで以下のことに気を付けてください。dtminlはhostで初期化されてdeviceでreductionされてhostに戻されています。このパート、別にhostと通信しなくても良いのですが、このパートの通信はどう書いても無くせませんでした。MPIminfind他のMPIprocessと通信するパートですが、hostで実行されるので特にGPU化を気にする必要はありません。hostのMPIで全体の最小のdtmingを求めてからdtを計算し、update toでdeviceに持っていくというアルゴリズムになっています。

hydro.cpp

void ControlTimestep(const GridArray<double>& G){
  using namespace mpi_config_mod;
  const double huge = 1.0e90;
  double dtminl = huge;
#pragma omp target teams distribute parallel for reduction(min:dtminl) collapse(3)
  for (int k=ks; k<=ke; k++)
    for (int j=js; j<=je; j++)
      for (int i=is; i<=ie; i++) {
	double ctot = sqrt(     P(ncsp,k,j,i)*P(ncsp,k,j,i)
			   + (  P(nbm1,k,j,i)*P(nbm1,k,j,i)
			      + P(nbm2,k,j,i)*P(nbm2,k,j,i)
			      + P(nbm3,k,j,i)*P(nbm3,k,j,i) ) / P(nden,k,j,i)); 
	double dtminloc = std::min({ (G.x1a(i+1)-G.x1a(i))/(std::abs(P(nve1,k,j,i))+ctot)
				    ,(G.x2a(j+1)-G.x2a(j))/(std::abs(P(nve2,k,j,i))+ctot)
				    ,(G.x3a(k+1)-G.x3a(k))/(std::abs(P(nve3,k,j,i))+ctot)});
	dtminl = std::min(dtminloc,dtminl);
	
      }
  //  Here dtminl is in host
  double dtming;
  int myid_wg;
  MPIminfind(dtminl,myid_w,dtming,myid_wg);
  dt = 0.05e0*dtming;
#pragma omp target update to (dt)
}

MPIによるゴーストメッシュの通信

CfCAではGPU4枚をMPIで通信しながら利用することができます。MPIの準備でちょっと特殊なのは以下の部分です。omp_get_num_devices()で使えるGPUの枚数が得られ、mpi_processごとにomp_set_default_device(gpuid)で違うGPUを割り当てます。

mpi_config.cpp

  // Select GPU device (OpenMP)
  ngpus = omp_get_num_devices();
  if (myid_w == 0) std::printf("num of GPUs = %d\n", ngpus);
  gpuid = (ngpus > 0) ? (myid_w % ngpus) : -1;
  if (gpuid >= 0) omp_set_default_device(gpuid);

流体計算ではゴーストメッシュ用のデータをパックしてMPI通信してアンパックするのが基本です。パックパートでは実メッシュの左端をBs.Xeという構造体につめて右端をBs.Xsに詰めます。その後MPI通信して左側ゴーストメッシュのBr.Xsと右側ゴーストメッシュのBr.Xeを受け取り、primitive variableのPに代入します。

boundary.cpp

void SetBoundaryCondition(FieldArray<double>& P,BoundaryArray<double>& Bs,BoundaryArray<double>& Br) {

  // x-direction
  // |     |Bs.Xe   Bs.Xs|     |
  // |Br.Xs|             |Br.Xe|
#pragma omp target teams distribute parallel for collapse(4)
  for(int n=0; n<nprim; n++)
  for(int k=ks; k<=ke; k++)
  for(int j=js; j<=je;j++)
  for(int i=0; i<ngh; i++) {
	  Bs.Xs(n,k,j,i) = P(n,k,j,ie-ngh+1+i);
	  Bs.Xe(n,k,j,i) = P(n,k,j,is      +i);
  }
 
  SendRecvBoundary(Bs,Br);

  
  // |     |Bs.Xe   Bs.Xs|     |
  // |Br.Xs|             |Br.Xe|
#pragma omp target teams distribute parallel for collapse(4)
  for (int n=0; n<nprim; n++)
  for (int k=ks; k<=ke; k++)
  for (int j=js; j<=je;j++)
  for (int i=0; i<ngh; i++) {
	  P(n,k,j,is-ngh+i) = Br.Xs(n,k,j,i);
	  P(n,k,j,ie+1  +i) = Br.Xe(n,k,j,i);
  }

};

MPIでの通信は以下のように行います。一般にMPI関数をhostで呼んで、送受信する配列のアドレスにGPU側のものを使います。通常は#pragma omp target use_device_addr(Bs.Xs_data)と書きたいところなのですが、この関数の引数にメンバー変数を書けないみたいなのでomp_get_mapped_ptr(Bs.Xs_data, dev);でdevice側のアドレスを持ってきています。ここでdevは最初に設定したGPU deviceのindexです。

boundary.cpp

void SendRecvBoundary(const BoundaryArray<double>& Bs,BoundaryArray<double>& Br){
  int dev = omp_get_default_device();
  //printf("myid gpuid=%i %i",myid_w,dev);
  int rc;
  int nreq = 0;
  
  if(ntiles[dir1] == 1){    
  // |     |Bs.Xe   Bs.Xs|     |
  // |Br.Xs|             |Br.Xe|
#pragma omp target teams distribute parallel for collapse(4)
    for (int n=0; n<nprim; n++)
    for (int k=ks; k<=ke; k++)
	for (int j=js; j<=je;j++)
	for (int i=0; i<ngh; i++) {
	    Br.Xs(n,k,j,i) = Bs.Xs(n,k,j,i);
	    Br.Xe(n,k,j,i) = Bs.Xe(n,k,j,i);
	  }
  } else {
  // |     |Bs.Xe   Bs.Xs|     |
  // |Br.Xs|             |Br.Xe|
    void* d_Bs_Xs = omp_get_mapped_ptr(Bs.Xs_data, dev);
    void* d_Bs_Xe = omp_get_mapped_ptr(Bs.Xe_data, dev);
    void* d_Br_Xs = omp_get_mapped_ptr(Br.Xs_data, dev);
    void* d_Br_Xe = omp_get_mapped_ptr(Br.Xe_data, dev);
    rc = MPI_Irecv(d_Br_Xs,Br.size1, MPI_DOUBLE, n1m, 1100, comm3d, &req[nreq++]);
    rc = MPI_Isend(d_Bs_Xe,Bs.size1, MPI_DOUBLE, n1m, 1200, comm3d, &req[nreq++]);
    rc = MPI_Irecv(d_Br_Xe,Br.size1, MPI_DOUBLE, n1p, 1200, comm3d, &req[nreq++]);
    rc = MPI_Isend(d_Bs_Xs,Bs.size1, MPI_DOUBLE, n1p, 1100, comm3d, &req[nreq++]);
   
   }

    if(nreq != 0) MPI_Waitall ( nreq, req, MPI_STATUSES_IGNORE);
    nreq = 0;	
};

まとめ

これまで多くの流体コードを書いてきたのですが、C++でGPU対応OpenMPはなかなか難しかったです。まずOpenMPがすぐhostとdeviceの通信をしたがるのでそれを抑制するのが大変した。関数の引数に主要な配列を書くこと、関数の中で大きな配列を確保しないことで対応しています。またC++でクラス、クラスのメンバー変数、メンバー関数を区別してGPU対応のためにはなるべく生の変数を書くというところ、omp target use_device_addrではメンバ変数が書けないところなど、やや面倒でした。とはいえ、ひとまとまりの情報を提供できたと一旦満足しています。