（２０２１年 9 月 2 日に開催されたアクセラレーション研究会における吉川氏（筑波大学）による応募書類の例です．参考にしてください．）--氏名：吉川 耕司
所属：筑波大学 計算科学研究センター課題名：ブラソフシミュレーションを用いた宇宙大規模構造形成におけるニュートリノの数値シミュレーショ
ン
1.

シミュレーションの科学的意義と手法
本応募において行った数値シミュレーションは，cold dark matter（以下、CDM）と質量をもつニュートリノの力学的影響を考慮した宇宙大規模構造の数値シミュレーションである。ニュートリノは宇宙年齢が 1 秒程度という非常に初期の段階で他の成分から脱結合した素粒子で，平均してニュートリノ 1 種類あたり 100 cm-3 という光子と同程度の高い数密度で宇宙全体に存在している。素粒子の標準模型ではニュートリノの質量は 0 とされていたため宇宙大規模構造の形成にはほとんど影響を及ぼさないと考えられていたが、スーパーカミオカンデなどによるニュートリノ振動の発見によってニュートリノが 0 ではない質量を持ち、少なくとも 2 種類のニュートリノは現在の宇宙で非相対論的になっていることが明らかとなった。これにより、将来の大規模銀河サーベイプロジェクトによる宇宙大規模構造の観測から素粒子物理学にとって重要なニュートリノ質量や質量階層を測定する可能性が指摘され、その為の理論モデルを詳細な数値シミュレーションによって構築する必要性が高まっている。本応募で行う数値シミュレーションはそのような必要性に応えるために必要不可欠なものとなっている。
従来の宇宙大規模構造形成の数値シミュレーションでは、N 体シミュレーションが事実上のデファクトスタンダードとして用いられてきたが、N 体シミュレーションは空間 3 次元と速度空間 3 次元からなる 6 次元
位相空間における物質分布をモンテカルロ的に質点（超粒子）でサンプリングするため数値シミュレーシ
ョン結果に必然的にショットノイズが内在する。特に、ニュートリノの様に質量が軽く速度分散が大きい成分は位相空間上での分布が広がっているため、CDM のように速度分散が極めて小さい成分と比較してショットノイズの影響がかなり大きくなってしまう。本応募で行った数値シミュレーションでは、このような N 体シミュレーションの欠点を克服する手法として、位相空間上の分布関数そのものをメッシュで離散化し、離散化された分布関数の時間発展を無衝突ボルツマン方程式（ブラソフ方程式）に従って数値シミュレーションする「ブラソフシミュレーション」という手法を採用する。
本数値シミュレーションでは、CDM と質量を持つニュートリノの 2 成分の重力相互作用のみを考慮し、宇宙論的な共動座標系で周期的境界条件の下でこの 2 成分の運動を数値的に解く。CDM は速度分散が極めて小さく位相空間での広がりが小さいため、N 体シミュレーションを用いてその運動を計算し、ニュートリノ成分については上で述べた理由によりブラソフシミュレーションの手法を採用する。したがって、
本数値シミュレーションは N 体シミュレーションとブラソフシミュレーションを組み合わせたハイブリッド計算を行っている（Yoshikawa et al. 2020）。
CDM の運動を計算する N 体シミュレーションについては、TreePM 法 （Bagla 2002） を用いた。

重力相互作用のうち周期境界条件下での長距離力については、Particle-Mesh 法を用いて重力ポテンシャルを Poisson 方程式に従って高速フーリエ変換を使って解くことによって求める。短距離力について
は、ある一定の半径内に存在する粒子からの重力を Tree 法で計算する。
ニュートリノの運動を計算するブラソフシミュレーションについては、6 次元位相空間の位置空間・速度空間それぞれを等間隔な 3 次元メッシュで離散化し、ブラソフ方程式を 6 つの次元について方向分割した 6 本の 1 次元移流方程式として有限体積法で時間発展させる。有限体積法のスキームは、我々が開発した空間 5 次精度の保存型セミラグランジュ法と単調性と正値性を維持する流束制限関数を組み合わせたものを用いた （Tanaka et al. 2017）。
2.

ハードウェア及びソフトウェア
本シミュレーションを行ったのは、理化学研究所 計算科学研究センターに設置されているスーパーコンピュータ「富岳」である。富岳に採用されているプロセッサは富士通製の A64FX で 48 個の演算コアを備えている。A64FX では Scalable Vector Extension （SVE） という SIMD 拡張命令を利用することができ、512 ビット幅の SIMD 演算用のレジスタを演算コアあたり 32 個備えている。したがって、16個の単精度演算、8 個の倍精度演算を同時に SIMD 実行することが可能である。A64FX の TDP は
約 200W である。
本シミュレーションにおける演算はその 95%以上がブラソフシミュレーションにおける移流方程式の数値計算と CDM の TreePM 法による N 体シミュレーションの粒子間相互作用の計算である。ブラソフシミュ
レーションにおける移流方程式の数値計算はそのほぼ全てを SVE 拡張命令セットによる単精度演算での SIMD 演算によって行っている。実際の実装には C 言語から SVE 命令を直接呼び出す Arm C
Language Extensions （ACLE） を用いた。また、TreePM 法による N 体シミュレーションの粒子間相互作用部分に関しても、x86 系のプロセッサで AVX-512 命令を用いて実装された PhantomGRAPE （Tanikawa et al. 2013）と呼ばれる粒子間相互作用の高速計算ライブラリを A64FX プロセッサ向けに SVE 拡張命令セットを用いて移植したものを利用した。
3.

電力性能
本シミュレーションの演算性能と電力性能を N 体シミュレーションとブラソフシミュレーションに分けて以下の様に求める。
TreePM 法による N 体シミュレーションでは、粒子間に働く重力を計算する演算が大部分を占めるのでこの演算性能を評価する。重力の計算では逆二乗則にカットオフ関数を掛けたものをテーブル参照で求めている。この計算には 1 相互作用あたり単精度浮動小数点で逆二乗則の計算に 38 演算とカットオフ関数の計算に積和演算 8 回（16 演算）が必要で、A64FX のコアあたり毎秒 1.2×109 相互作用の計算を達成している。従って、N 体シミュレーションの演算数は単精度演算でプロセッサあたり(1.2 × 109 ) × (38 + 16) × 48 = 1.56 × 1012

回で、演算性能は 1.56Tflops である。また、当該部分でのプロセッサあたりの電力消費量は富岳で利用できる富士通製の CPU 性能プロファイラの測定結果によると 144W であり、電力性能は単精度演算であることによる補正を行うと 5.42Gflops/W である。
一方、ブラソフシミュレーションについては、位相空間の各次元に沿っての移流方程式の数値解法がほぼ 100%を占める。富岳で利用できる富士通製のプロファイラを用いた測定から、該当箇所のプロセッサあたりの演算性能と電力消費量を運動量空間と座標空間での計算に分けて求めると運動量空間では演算性能は 920Gflops で電力消費量は 110W、座標空間では演算性能は 600Gflops で電力消費 量 は 120W で ある 。 従って 、単 精 度 演 算 であ る ことによる 補 正 後 の 電力 性能 はそ れ ぞれ 、
4.18Gflops/W と 2.50Gflops/W である。
シミュレーション全体における N 体シミュレーションとブラソフシミュレーションの運動量空間と座標空間の実行時間の割合はそれぞれ 25%、35%、40%であるので、シミュレーション全体での電力性能は実行時間で重み付けをして 3.81Gflops/W とした。
参考文献
Bagla, J.S. 2002, Journal of Astrophysics and Astronomy, 23, 185
Yoshikawa, K., Tanaka, S., Yoshida, N., Saito, S. 2020, ApJ, 904, 159
Tanaka, S., Yoshikawa, K., Minoshima, T., & Yoshida, N. 2017, ApJ, 849, 76
Tanikawa, A., Yoshikawa, K., Nitadori, K., Okamoto, T. 2013, New Astronomy, 19, 74