導入
HPCの世界情勢としてGPUを使うことが求められているのですが、現状、Fortran+OpenMPだとうまくいかないことがわかりました。OpenACCが使える環境であれば良いのですが、NVIDIAの環境に依存してしまいます。そこでメーカーに依存しないコードを求めて、FortranのコードをC++に移植してみたいと思います。
長いコードの移植は大変で現実的ではないとも思っていたのですが、一つ予想外だったのはAIが書いたコードを使えることが多く、AIに相談しながらやれば移植も無理ではないかもしれないと思えてきました。
FortranユーザーがまずひっかかるのはCでの多次元配列の扱いです。今回やってみて、やっぱり難しいなと思いましたが（後述）、数年前よりは状況がいいなと思いました。というのも今回chatgpt(有料版)に「Athena++を参考にして多次元配列を扱うC++のヘッダファイルを書いてください。」と命令したらちゃんと動くコードを書いてくれました。実際には、GPU側の都合でそうすんなりとはいかず、苦戦する部分もあったのですが、AIは強い味方ではあります。
実装の解説
先ほどヘッダをAIに書かせたと書きましたが、流体計算のコードそのものもAIが書いてくれました。少なくともコンパイルは通りました。今回紹介するものはそこから修正したものです。
全体の構造
全体の構造はこんな感じです。昔書いたFortran + OpenACCのコードと似たような構成となっています。
int main() {
using namespace resolution_mod;
using namespace hydflux_mod;
using namespace boundary_mod;
printf("setup grids and fields\n");
AllocateHydroVariables(U,Fx,Fy,Fz,P);
AllocateBoundaryVariables(Xs,Xe,Ys,Ye,Zs,Ze);
printf("grid size for x y z = %i %i %i\n",nx,ny,nz);
GenerateProblem(P,U);
printf("entering main loop\n");
int step = 0;
auto time_begin = std::chrono::high_resolution_clock::now();
for (step=0;step<stepmax;step++){
ControlTimestep();
SetBoundaryCondition(P,Xs,Xe,Ys,Ye,Zs,Ze);
GetNumericalFlux1(P,Fx);
GetNumericalFlux2(P,Fy);
GetNumericalFlux3(P,Fz);
UpdateConservU(Fx,Fy,Fz,U);
UpdatePrimitvP(U,P);
time_sim += dt;
//printf("dt=%e\n",dt);
if (step % stepsnap == 0) {
#pragma omp target update from (P.data[0:P.size])
Output1D();
}
}
auto time_end = std::chrono::high_resolution_clock::now();
std::chrono::duration<double> elapsed = time_end - time_begin;
printf("sim time [s]: %e\n", elapsed.count());
printf("time/count/cell : %e\n", elapsed.count()/(nx*ny*nz)/stepmax);
printf("program has been finished\n");
return 0;
}
今回は一番簡単な流体計算として１次精度の風上法を実装してます。精度をあげたり、数値フラックスをHLLEで計算したものについてはこちらを参照のこと。実装によっては空間のindexが入った中間変数が必要になりますが、なるべくそうはしないほうが良いです。コンパイラがアホな時、空間のindexが入った中間変数をホストと通信しようとします。
void GetNumericalFlux1(const Array4D<double>& P,Array4D<double>& Fx){
#pragma omp target teams distribute parallel for collapse(3)
for (int k=ks; k<=ke; ++k)
for (int j=js; j<=je; ++j)
for (int i=is; i<=ie+1; ++i) {
double qL = P(nden,k,j,i-1);
double qR = P(nden,k,j,i );
double ux = 0.5e0*(P(nvex,k,j,i-1)+P(nvex,k,j,i));
double q_up = (ux >= 0.0) ? qL : qR;
Fx(mden,k,j,i) = ux * q_up;
}
}
ここでP, Fxを関数の引数にしています。Fortranユーザーからすると参照渡し以外を使うのは混乱の元なのでArray4D<double>&のように常に&をつけましょう。この関数内でPは変更しないのでconstがついてます。
P, Fxは実はnamespace hydflux_modで定義しており、関数の引数にする必要は本来ないのですが、引数にしてあげたほうがコンパイラが全体の構造を把握しやすいみたいです。引数にしないときにはこのループの前でホストからデバイスにデータを送信して、後にデバイスからホストにデータを送信していましたが、このようにすると無駄なホストとの通信がなくなりました。ちなみにnamespaceはFortranで言うところのmoduleです, hydflux_modのmodはそこからきています。
GPU化に関してはOpenMP 4の#pragma omp target teams distribute parallel forを使用しています。OpenMP 5では#pragma omp loop order(concurrent)が使えるわけですが、むしろ
最適化されない場合もあるので、わざわざ古い書き方をしています。
データの取り扱いの一般論
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])
多次元配列の取り扱い
先ほどのhydro.cppでP(n,k,j,i)はprimitive 変数の配列（的なもの）です。P(nden,k,j,i)は密度です。P(nvex,k,j,i)はx方向速度です。この変数はiのところが連続になるようメモリを確保しております。こういう確保の仕方をSoA(Structues of Array)といいます。逆に密度、速度などが連続になるような扱いをAoS (Array of Structures)と呼びます。どちらが良いのかは場合によるかと思いますが、一般にSoAが有利と言われているようです。
C言語の多次元配列は本来P[n][k][j][i]みたいな感じですが、()のオペレータを定義して、直接アクセスできるようにしています（上で配列的なものと書いたのはそういう意図）。Pそのものは構造体で１次元化されたデータと自分の次元の情報を持っています。このパートを自分でで書くのは結構難しく（例えばこの記事をご覧ください。）、まぁ呪文だと思って他の人のものをコピーしたり、AIに書いてもらうと良いです。
この多次元配列の実装はもっと高尚にもできるのですが、GPU対応という観点からはむしろしないほうが良いです。例えばdataとかデータの数などをprivate変数にして見えなくすることもできますが、GPUに乗せるにあたって必要なこともあるのでわざわざ見えるようにpublicにしています。またstd::vectorを使うこともできるのですが、バグった時に何をやってるのかわかりにくくなると思いました。このコードのOpenACC版では使用しているので興味がある方はご覧ください。
template <typename T>
class Array4D {
public:
T* data = nullptr;
int nv = 0, n3 = 0, n2 = 0, n1 = 0;
size_t size = 0;
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];
}
};
多次元配列のallocateは以下のようにしています。GPUに乗せる時に難しさがあります。このUは構造体であり、配列ではありません。構造体のメンバーの変数(U.data, U.n1など)は配列だったりスカラーだったりしますのでOpenMPの指示文に書くことができます（逆に言うとここに構造体のUを書いてもdevice側にどのていどallocateされるのか不明）。手順にそってGPU側でallocしてからtoでmacしています。ここを飛ばしてallocせずにmapしたら上手くいかなかったです。気をつけてください。このように構造体を作ってもompの指示文に書く時には生の配列が必要になります。
void AllocateHydroVariables(Array4D<double>& U){
U.allocate(mconsv,ktot,jtot,itot);
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 enter data map (alloc: U.data[0: U.size])
#pragma omp target update to ( U.data[0: U.size], U.n1, U.n2, U.n3, U.nv)
}
ヘッダファイルの書き方
C++では宣言と定義を分けて考え、宣言はhppファイルに定義はcppファイルに書きます。Fortranユーザーからすると「なんで同じようなことを２回書く？」と思えるが、これは言語の仕様の重要な部分だそうです。
宣言は他のファイルからアクセスするためのサマリーのように考えるとよく、具体的には宣言にはexternを用います。
コードの書き方として解像度以外の部分はhppファイルで宣言、cppファイルで定義という書き方を保っています。宣言の場合はexternを用います。
#ifndef HYDRO_HPP_
#define HYDRO_HPP_
#include <cmath>
#include <cstdio>
#include <algorithm>
#include "hydro_arrays.hpp"
namespace hydflux_mod {
extern int mconsv; //!
extern Array4D<double> U; //! U(mconsv,ktot,jtot,itot)
extern int mden,mrvx,mrvy,mrvz,meto;
extern Array4D<double> Fx,Fy,Fz;
extern int nprim; //!
extern Array4D<double> P; //! P(nprim,ktot,jtot,itot)
extern int nden,nvex,nvey,nvez,nene;
};
void AllocateHydroVariables(Array4D<double>& U, Array4D<double>& Fx,Array4D<double>& Fy,Array4D<double>& Fz,Array4D<double>& P);
void DeallocateHydroVariables(Array4D<double>& U, Array4D<double>& Fx,Array4D<double>& Fy,Array4D<double>& Fz,Array4D<double>& P);
void GetNumericalFlux1(const Array4D<double>& P,Array4D<double>& Fx);
void GetNumericalFlux2(const Array4D<double>& P,Array4D<double>& Fy);
void GetNumericalFlux3(const Array4D<double>& P,Array4D<double>& Fz);
void UpdateConservU(const Array4D<double>& Fx,const Array4D<double>& Fy,const Array4D<double>& Fz,Array4D<double>& U);
void UpdatePrimitvP(const Array4D<double>& U,Array4D<double>& P);
void ControlTimestep();
#endif
ヘッダファイルは２回読み込まれないためのおまじないとして#ifndef HYDRO_HPP_みたいなものをつけます。
cppファイルでは実際の定義を行っています。
namespace hydflux_mod {
#pragma omp declare target
int mconsv{5}; //!
HydroArrays<double> U; //! U(mconsv,ktot,jtot,itot)
int mden{0},mrvx{1},mrvy{2},mrvz{3},meto{4};
HydroArrays<double> fluxx,fluxy,fluxz;
int nprim{5}; //!
HydroArrays<double> P; //! P(nprim,ktot,jtot,itot)
int nden{0},nvex{1},nvey{2},nvez{3},nene{4};
#pragma omp end declare target
};
解像度などの定数の取り扱い
resolution_modで解像度を設定しています。これはヘッダファイルで設定されています。先述の通り、C++の仕様としてはヘッダファイルは宣言だけして定義しないのですが、２回書くのが面倒ですし、cppファイルで定義すると複数のファイルから参照するが難しいので、ここで定義もしてしまっています。変数にはinlineを定数にはinline constexprを使います。
namespace resolution_mod {
inline constexpr int stepmax{100}; // max step
inline constexpr int stepsnap = stepmax/100;
inline double time = 0.0e0;
inline double dt;
inline constexpr double timemax = 5.0e0;
inline constexpr double dtout = 5.0e0/600;
inline constexpr int nx{64}; //! resolution for x
...
};
まとめ
GPUを使う時Fortranのコンパイラに不満がでがちです。大変ですがC++に移植してしまうのも一つの手段です。今回はFortranユーザーの気持ちに沿ってC++のコードの書き方を解説してみました。自分がつまづいたところを中心に説明してあります。全体の印象としてOpenACCのほうがやっぱり簡単でOpenMPは気が利かず発展途上のように感じます。コンパイラが成熟する日を願っています。