Skip to content

Instantly share code, notes, and snippets.

@kumpeishiraishi

kumpeishiraishi/speeding_up_your_particle_simulations.org

Last active October 25, 2023 19:26
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save kumpeishiraishi/1a0f3932adbc4044a152699c07dd4fc9 to your computer and use it in GitHub Desktop.
Save kumpeishiraishi/1a0f3932adbc4044a152699c07dd4fc9 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
speeding_up_your_particle_simulations.org

粒子系シミュレーションにおける高速化のメモ

個人的なチェックリスト

  1. まずは普通にコードを完成させる
    • MDなら、NVEで時間発展させて全エネルギーの保存を確認し、コードにバグがないことの確信を得る
      • 全エネルギーの揺らぎが、運動エネルギー・ポテンシャルエネルギーの揺らぎと比べて、2桁程度以上小さければOK
    • MCなら、MSDなどを計算して、自分のコードのダイナミクスと先行研究に齟齬がないことを確認する（その他の方法だと、サンプルした配置のポテンシャルの比較？）
    • この段階では、高速化の諸々を先取りするといった余計なことはしない
    • まずは、とにかく素直に正しいコードを書き終える方が、作業は結局速く終わる（経験上）
  2. アルゴリズム上の高速化手法を実装する
    • Allen-TildesleyやFrenkel-Smitといった教科書に載っている、ベルレの隣接リストやリンクセルのこと
    • デバッグのため、適当な長さのシミュレーションの後に終配置を出力する
    • この際、最初に書いた単純なコードと、出力される粒子位置が 完全に 一致することを確認（ diff などでチェックできる）
      • MCは（乱数のseedが同じなら） いつでも必ず 一致するはず
      • MDは、足し算の順番によって結果が微妙に変わることがあるが、参照する粒子番号の適切な並べ替えによって一致させることは可能
    • 経験上、粒子数が少ないとき（コード開発時はそうしていると思う）は、素朴に隣接リストを構築するのが速い（3次元のLennard-Jones相互作用で1000粒子程度）
    • 経験上、以下のケースでは、隣接リストの構築時にリンクセルを併用して相互作用粒子を探すのが速い
    • リンクセル+隣接リストのコードを書く際は、素朴な $ij$ ループで隣接リストを構築するステップを明示的に挟むと、デバッグしやすいかも
  3. ソフトウェアパイプライニングを施す
    • 渡辺先生の記事を参照
    • MDの場合は力計算カーネル、MCの場合は粒子 $i$ のポテンシャル計算カーネルについて、ループ順序を適切に変更する
    • 終配置の位置は 完全に 一致させる
  4. ループアンロール
    • ここまではSIMDのことは一切考えない
    • 以上でコード中のループ構造が確定するはずなので、ここでようやくSIMDについて考え始める
    • 最速のコードにならないままSIMDに関して試行錯誤し始めると、結局時間を無駄にすることになる
    • SIMD化の際はいきなりintrinsicを書き始めず、ループアンロールから始める（渡辺先生の記事を参照）
    • MDの場合は、 $j$ 粒子のループ（力計算カーネルの内側ループ）を4倍アンロール
    • MCの場合は、ポテンシャル計算カーネルのループを4倍アンロール（MCの場合、 $j$ 粒子のループしか回らないはずなので自明ですが）
    • ここでも終配置の位置は 完全に 一致させる
  5. SIMD化
    • 渡辺先生の記事を参照
    • MCの場合、終配置の位置は 完全に 一致させる
    • MDでも完全に一致させたい所だが、コンパイラによってはSIMDを使ったコードと計算結果が微妙に異なることがある
    • その場合は、元のコードとSIMD化コードで同一配置を読み込み、計算した力を（各粒子各成分毎に）書き出す
    • 各成分毎に両データの差を取ったとき、その最大値が $10-14$ 程度になることがある
    • 同様の理由で、コンパイラによる最適化に注意
      • SIMD化を一切考えない前ステップまでで、CPUタイプに適したコードを生成するオプション（GCCの場合 -march=native など）を付した場合、コンパイラによって自動でSIMD化されることがある
      • すると、同様の理由でMDの終配置が一致しない
      • 実際上は、SIMDコードを書く前はこうしたオプションを付与せずに終配置の完全一致を確認、SIMDコードを書く段になってオプション付与・力の一致を確認、としている
    • 同じこと（ソフトウェアパイプライニング、ループアンロール、SIMD化）を隣接リスト周りでもやると、もうちょっと速くなる

バグに遭遇したら

  • まず時間発展を止めて、単一配置の力・ポテンシャルの計算に限定する
  • 上で概観したのは力・ポテンシャル計算カーネルの書き換えであるが、問題は以下の2点に切り分けられるはず
    1. 相互作用する粒子対は正しく選択できているか？
    2. 選択した粒子対に関して、ポテンシャルは正しく計算できているか？
  • まず1点目をチェックするため、粒子0と相互作用する粒子の番号、粒子1と相互作用する粒子の番号、、、というファイルを出力し、元のコードと比較する（ diff でチェック）
  • ループが正しい粒子の組み合わせについて回っていることを確認したら、粒子間距離 $rij$ 、半径の和 $σij$ 、ポテンシャルや一階微分を出力し、元のコードと比較する

ご意見・ご感想はお気軽にお寄せください（https://kumpeishiraishi.github.io）

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment