SPH（粒子法）シミュレーションコードの GPU 化に取り掛かりました。
まずは、土のみ3Dコードから。テストには 斜面崩壊（粒子数77万弱、200,00step）を使用しました。実行環境は Dtransu と同じです。
GPU: NVIDIA RTX 4000 Ada(20GB、CC 8.9)
CPU: AMD EPYC 9754(128 コア)
いくつかの試行を得て、こちらは GPU 版が CPU 版に勝てない結果となりました。
OpenMP 32 / 64 / 128 スレッド 94.4 / 60.3 / 48.0s
OpenACC 75〜84s
CUDA Fortran + Thrust（set_box のみ) 41.7〜55.3s
計算コードの変数を倍精度で計算していたため、FP64律速かつメモリ帯域律速となり、搭載GPU にとって二重に不利、CPU に有利でした。特に処理時間の約45%を占める make_interaction とkernel_correction は、GPU 側で改善できませんでした。粒子のセル順ソートは GPU では悪化、コアレス化は微妙でした。
セル分割処理だけはGPUで効果がありました。この処理は整数演算 が主体で、FP64律速にも帯域律速にも縛られなかったのでしょう。
- セル分割（set_box）：空間を格子(セル)に区切り、「どの粒子がどのセルにいるか」の対応表（名簿）を作る処理、近傍探索の前準備です。GPU 化では内部のソートに Thrust の sort_by_key （CUB の基数ソート） を使い、インデックス（cell_id と粒子番号）だけを並べ替えて名簿を作ります（粒子データ本体は動かさない）。これが GPU で効いた処理です。
- 粒子のセル順ソート（reorder_particles）：粒子データ本体そのものを、セルの順番に並べ替えておく「席替え」(50 step 毎)。近くにある粒子をメモリ上でも近くに置き、キャッシュ効率を上げる狙い。set_box内のソートと違い、倍精度の粒子配列を丸ごと移動するため、GPU では悪化しました （CPU ではキャッシュ効果で改善）。
- メモリアクセスのコアレス化（pair_layout）：粒子ペアのデータを GPU が読みやすい並び順に変える。GPU は隣り合うスレッドが連続したメモリを読むと速い(=コアレス化)ので、その並びに合わせる。
セル分割＋ソート＋コアレス化は DualSPHysics の講義で教えていただいたのを覚えていました。当時はメモリ上の最適配置程度で大きな影響があるとは思わず、ちょっとしたテクニック程度の認識でした。ほかにも、DualSPHysicsでは dx/dw など過去の値に依存しない変数は FP32、累積する変数は FP64 でドリフトを避ける工夫がなされています。
土のみコードは両者を倍精度で計算しており、一部を単精度にすると過去の結果（降伏、流体化の判定）の完全再現が困難。単精度化は周囲の応力均衡を変化させ、隣接粒子の降伏タイミングを次々とずらしていく連鎖を引き起こします。これは避けたいので、今回は速度を犠牲にして見送り。
今回はここまで。次は、残したFP32化に手を付けるか、良いGPUに入れ替えるか、でしょうか。A100 / H100 クラスの GPU であれば、CPUを超えるかもしれません。
