特集続スーパーコンピュータ 京 の利用 図 -1 シミュレーションによって再現された, ダークマター構造の進化の様子. 明るさはダークマターの空間密度を表し, 明るいところは密度が高い. 宇宙初期 ( 一番左のパネル ) に存在した微小な密度揺らぎが, 時間経過とともに ( 左から右 ) 重力により

Transcription

1 特集続スーパーコンピュータ 京 の利用 基応専般 8 京 の威力で宇宙の正体に迫る ダークマターの超大規模シミュレーション 石山智明 筑波大学計算科学研究センター ダークマターとは 夜空に煌めく無数の星々や, 我々が住む地球の親星である太陽は銀河系 ( 天の川銀河 ) のメンバである. 銀河系は半径数万光年程度のディスク状に分布した数千億個の星々からなると考えられている. 銀河が多数群集まっているのが銀河団であり, それらはさらに超銀河団という集団構造を形成している. 超銀河団同士もフィラメントやシート状に分布する銀河でつながっていて, 宇宙の大規模構造を形成する. このように宇宙は階層的な構造をしている. こういった我々の目で観測できる天体は, 原子や分子のような バリオン と呼ばれる物質で構成されている. 最新の宇宙物理学は, バリオン物質は宇宙の構成要素全体の 4 ~ 5% 程度の割合でしか存在しないことを明らかにした. 残り 25% 程度は ダークマター と呼ばれる重力のみ作用する物質, 残りは宇宙の加速膨張にかかわる ダークエネルギー である. 宇宙誕生からおよそ 140 億年間にわたって, どのようにダークマターが重力的な構造形成の主要な役割を果たしてきたかを俯瞰してみよう. 宇宙誕生直後は, ダークマターは空間全体にほとんど一様に存在していた. しかしいたるところに, 密度が少しだけ高いところ, 低いところ ( 密度揺らぎ ) が存在した. ある程度高密度な領域では, 重力によって宇宙膨張に打ち勝つことができ収縮していく. こうしてできる安定なダークマター天体が, ダークマターハローである. そしてハローの中心の深い重力ポテンシャルの中にバリオンが集まる. バリオンが高密度に凝縮することで原始星ができ, ハローが合体を繰り返すことで, 現在観測される銀河のような多種多 様な天体を形成していったと考えられている 1. 図 -1 にシミュレーションによって再現された, 宇宙のダークマター大規模構造が形成し進化する様子を示した. ダークマターの存在は, さまざまな観測で間接的に確かめられている. 間接的というのは, 星や銀河のように光を介して直接観測できないが, 見える物質 による重力だけでは説明できない観測結果が数多く存在するという意味である.1 つの例として, 回転曲線問題というものがある. 銀河の星々は銀河中心に対しておおよそ回転運動をしている. 回転速度はその星の位置より内側にある物質の総質量によっていて, 物質が多いほど速度も大きい. ところが, 観測された星の回転速度は, 星の総質量から推定されるものより大きい. これは見えない物質, つまりダークマターが存在することを強く示唆している. 逆に回転速度からダークマターの量を見積もることができ, 銀河系は半径数十万光年にわたって広がった,1 兆太陽質量程度の, 楕円状のハローが作る巨大な重力場の中心に存在していると考えられている. 理論的にも, バリオンの密度ゆらぎの大きさだけでは, 現在の宇宙の構造を作るには不十分であることが分かっている. よってダークマターは天体形成に必要不可欠であり, 逆にダークマターの分布がバリオンの集まり方, つまり銀河や銀河団などの分布を決定する. その分布には宇宙初期の密度揺らぎの性質についてや, 宇宙年齢, 各物質成分の割合などの情報が刻み込まれている 1). それでは宇宙の構造形成, 進化に重要な役割を担ってきたダークマターとは一体何だろうか? 実は 1 観測される銀河は, バリオン質量の数倍から数十倍にもおよぶダークマターハローの中心に存在していると考えられている. 829

2 特集続スーパーコンピュータ 京 の利用 図 -1 シミュレーションによって再現された, ダークマター構造の進化の様子. 明るさはダークマターの空間密度を表し, 明るいところは密度が高い. 宇宙初期 ( 一番左のパネル ) に存在した微小な密度揺らぎが, 時間経過とともに ( 左から右 ) 重力により発達していく. 一番右のパネルは宇宙誕生から約 140 億年後の現在の宇宙の姿を表している. 図中右下の数値は空間スケールを表し,1Mpc はおよそ 300 万光年である. いたるところに存在する楕円状の高密度領域それぞれが個々のダークマターハローである. ダークマターの素粒子としての正体はよく分かっておらず, 宇宙物理学的にも, 素粒子物理学的にも重要な未解決問題の 1 つである. 有力な候補として 超対称性粒子 である ニュートラリーノ が提案されている. ニュートラリーノは, 自己対消滅をしガンマ線を放出するという特徴がある. このガンマ線を検出することで, 間接的にダークマターを検出しようとする試みが世界中でなされているが, 決定的な結果は得られていない. スーパーコンピュータ シミュレーション 宇宙の構造形成過程は, 重力の非線形性が本質的に重要な役割を果たしているため, そのダイナミクスを研究するには数値シミュレーションが非常に有用である. 宇宙空間の物理現象は, 地上で実験して検証することが不可能であり, また質量, 空間, 時間スケールが非常に広大である. したがって, スーパーコンピュータによる大規模シミュレーションが必要不可欠である. 時にシミュレーションは 理論の望遠鏡 とも呼ばれ, 宇宙の研究において非常に重要な位置を占めている. 特にダークマターは直接観測できないので, ダークマター構造形成の研究は, シミュレーションによって大きく進められてきた. 手法としては, 重力多体シミュレーション, いわゆる N 体シミュレーシ ョンがよく用いられる. これは対象を N 個の粒子としてモデル化し, 粒子間にはたらく相互作用重力を計算し, 運動方程式を解いて, 粒子の時間発展を追うシミュレーション手法である. ダークマターのほかに, 星団や銀河, ブラックホールなど適用範囲は非常に幅広い. ダークマターの場合, 個々の粒子は位相空間上のダークマターを離散化したものをあらわす. シミュレーション粒子がダークマター素粒子に相当するわけではない. 宇宙初期の密度揺らぎを持ったダークマターの分布からスタートし, 重力相互作用による時間発展を追うものを, 特に宇宙論的シミュレーションと呼ぶ. 無限遠方からの重力を取り入れるために周期境界条件が用いられる. 初期密度揺らぎは宇宙背景放射の観測などから非常によく制限されている. 対象とする天体によって, 扱う領域はさまざまである. 銀河の空間分布を調べたい場合は, 観測できる宇宙の全領域にせまる領域を扱う. 銀河の微細構造を調べたり, 我々の銀河の形成過程を研究するときは 1 つのハローだけを扱うことが多い. シミュレーションに用いる粒子数はきわめて重要である. どれだけ広い宇宙をシミュレーションできるかが決まり, 銀河分布の精密さに影響する. またどれだけ微細なダークマター構造が分解できるかが決まり, それは銀河系自身のモデルや, ダークマター対消滅シグナルの予測の正確さに直結する. 830

3 Tree 図 -2 TreePM 法の概念図. 粒子分布を 8 分木で表現する ( 左パネル ). ツリー法で計算するのは, ある粒子から見てカットオフ半径 ( 緑色の円 ) 内にある粒子からの力のみである. それより外側は PM のみで計算される. ツリー法では, 十分に離れた粒子群からの力はまとめて多重極展開として計算される ( 緑丸 ). 右下のグラフは粒子間距離の関数としての, 各力の大きさである. 緑線がツリー法で, 青線が PM で計算される力である. これらの合計がニュートン重力 ( 破線 ) である. 重力はすべての粒子ペア間にはたらくが, 粒子数 N に対して, ペア数は N 2 に比例する. そのまま計算すると, 粒子数を 10 倍にしたときに計算量は 100 倍となる. これでは粒子数を増やすのは非常に困難である. そこでツリー法 2) のような近似法がよく使われる. 基本的なアイディアは, 近傍の粒子からの力は直接計算するが, 十分遠方の粒子群からの力はまとめて多重極展開からの力として計算するというものである. ツリー法では計算量を O (N log N) に減らすことができ, 粒子数を大きく増やすことが可能となる. 宇宙論的シミュレーションではツリー法をいわゆる Particle-Mesh(PM) 法と組み合わせた TreePM 法がよく用いられる ( 図 -2). この方法では, ツリー法で計算するのは, ある粒子から見てカット オフ半径内 2 にある粒子からの力のみである. カットオフ半径より外側は,PM 力からの寄与のみとなる. カットオフの導入により, 通信量を大きく減らすことができる.PM 法では, まず粒子の位置と質量から一様メッシュ状の密度場を計算し,FFT とポアソン方程式から一様メッシュ状のポテンシャルを計算する. そして差分と補間から, 各粒子位置 への力を計算する. また周期境界条件を自然に組み込むことができる. 宇宙論的シミュレーションでは, 前述のように大規模化が必要不可欠であるため,TreePM 法を用いた大規模並列コードの開発が世界のいくつかのグループで進められてきた. しかし, これまで約 1,000 並列度程度までに対応した並列化が報告されているくらいであった. 我々のグループでも C++ による OpenMP+MPI 並列コード,"GreeM" を開発してきた. 京 のような数千並列以上の大規模環境での代表的なボトルネックには, (1) 構造の時間発展によるロードインバランス (2)PM 法における, メッシュ構造の大規模全対全 通信 3 などがあるが, 我々は新しい並列化手法を開発し,8 万ノードにもおよぶ 京 のフルシステム上で効率良く動作するアプリケーションの開発に成功した. 開発にあたっては,2014 年 2 月現在, 理化学研究所計算科学研究機構粒子系シミュレータ研究チームの似鳥啓吾研究員, 同チームリーダーの牧野淳一郎教授の協力を得た. 構造の発達によるロードインバランス並列化は, シミュレーション空間を並列数で分割し, 各ノードに割り当てることで行う. 宇宙論的シミュレーションでは図 -1 から分かるように, 時間発展とともにいたるところに高密度なダークマターハローが発達する. よって領域分割が空間一様で時間不変である場合, 構造が発達するとロードバランスが非常に悪化する. 多くのコードでは,morton ordering や peano hirbert curve のような空間充填曲線に沿って粒子を並べ, 並列数で均等に分割して領域分割を行っている. この方法ではロードバランスを均等にしやすい反面, 領域形状が複雑化し, ステップ間での粒子交換や, ツリー法による力の計算に必要な, 他のノードに存在するツリー構造の通信粒度が細かくなる傾向にある. 8 京 の威力で宇宙の正体に迫る ダークマターの超大規模シミュレーション 2 許容できる誤差と計算量の観点から, 典型的には PM で用いるメッシュ幅の 3 倍程度に設定する. 3 通信にかかわる全ノードから全ノードへのデータの相互交換. 831

4 特集続スーパーコンピュータ 京 の利用 これらの方法により, 各ステップで動的に, ロードバランスを調整しながら領域分割を行っている. さらにサンプリングによる浮動を最小限にするために, 過去数ステップの領域形状を平均化している. こうして, 最小限の通信で 8 万ノードでも計算負荷をほぼ均等に保つことに成功した. 図 -3 領域分割の例.2 次元で 8 8 分割している. 右上のパネルは中心部を拡大したものであり, 高密度なハロー部分で細かく分割されていることが分かる. そこで我々のコードでは, 再帰的多段分割法に基づいた 3 次元領域分割を採用した.3 軸それぞれの分割数を決め, まず x 軸方向に計算時間が均等になるように分割する. そして分割された空間ごとに独立に,y 軸方向,z 軸方向に分割していく. 結果領域形状は直方体となり, 通信粒度が小さくなることを防ぐことができる. 領域分割の例を図 -3 に示す. 構造が発達する個所で細分化されている様子が分かる. またこの領域分割は 京 のようなトーラスネットワークにマップしやすく, 通信最適化を促進しやすいというメリットがある. 通信量を抑えるために, 領域分割はサンプル粒子に対して行っている. 各プロセスは適当なサンプリングレートで粒子をサンプルし, ルートノードがそれらを集める. そして領域再分割を行い, 結果を全体に送りかえす. このとき, 各ノードは 1 つ前のステップの計算時間に比例するように, サンプリングレートを決める. 前のステップで少し計算負荷が大きかったノードは, サンプルされる粒子数も少し大きくなる. 領域再設定は, 各ノードでのサンプル粒子の数が均等になるように行われる. したがってサンプル粒子数が大きかったノードは, 次のステップで計算負荷が減るように領域再設定が行われる. FFT 計算のための, メッシュ構造の大規模全対全通信 PM で力を計算するとき FFT を実行するが, 我々のコードでは 1 次元スラブ並列にのみ対応した FFTW ライブラリ 4, または 3 次元並列に対応した富士通製の FFT ライブラリを使用している. いずれの場合でも,FFT 用の分割されたメッシュデータの構造は, ノード間で均一でなければならない. ところが領域形状は直方体であり, 図 -3 でも見られるように,3 軸不等の度合いがかなり大きくなることもある. したがって各ノードの担当領域を覆うローカルなメッシュは, ノード間でサイズや形状が大きく異なるが, それを均一なスラブ分割, または 3 次元分割になるようにデータ構造を変換しなければならない (FFT で計算されたメッシュ上のポテンシャルに対しても同じ変換が必要 ). 計算全体の中で FFT 自体の計算量はたいしたことはないが, このデータ構造を変換する通信は大規模で, 非常に不均一な全対全通信であるため, 数万並列規模だとボトルネックになり得る. 我々は全体をいくつかのグループに分け, まずグループ内での全対全通信を行ってからグループ間で通信するというように,2 段階に分けて通信を行うことによって性能を大きく改善することができた. 基本的に FFT を実行するノード数は, 全体のノード数より小さい ( 特に 1 次元スラブ並列にのみ対応した FFTW を用いる場合 ). まず, 各グループは FFT を実行するノード全体を内包できるサイズ, 形状にし,FFT 用のメッシュ全体をグループ内で分散して持つようにする. そして, グループのメンバノードが持つ部分メッシュを, グループ内全

5 対全通信で FFT 用のデータ構造へと変換する. この時点では, 各グループは全体の FFT 用メッシュを持っているものの, 中身はグループが覆う領域からの寄与のみなので不完全である. 次にグループ間通信を行い, 各メッシュ上の値を足し合わせることで, 全体の完全な FFT 用メッシュを生成する. こうすることで, この通信部分を数倍高速化することに成功した 3). これらの並列化手法を実装することで, 京 のフルシステム上でもほぼ完璧なロードバランスを実現した. さらに粒子間の相互作用重力を高速に計算することが重要である. 富士通コンパイラ提供の組込み関数を用いて手動 SIMD 化を行い, さらに力を受ける粒子とおよぼす粒子両方のループで, 手動ループアンロールを行った. 最適な展開数を見つけるために, 展開数を指定すれば自動でコードが生成されるようにした. 我々の実装では,1 相互作用の計算には 51 演算必要であり,FMA 部と非 FMA 部の割合が 1:1(17 命令ずつ ) であるので, 理論ピーク性能は CPU ピーク性能の 75% となる. この部分だけの単純なベンチマークテストでは,72% 程度の実行効率 ( 理論ピークの 96%) を達成した ( この部分は似鳥氏の実装による ). 図 -4 に GreeM のストロングスケーリングの測定結果を示した. 粒子数 1 兆のシミュレーションでは, 京 のフルシステムまでほぼ理想的にスケールしていることが分かる. 最終的な実効性能は,2 兆ダークマター粒子のシミュレーションを, 京 のほぼ全システムを用いて,5.67 ペタフロップスの実効性能 ( 実行効率 55%) で達成した. この成果が認められ, ハイ パフォーマンス コンピューティングに関する国際会議 SC 12 ( 2012 年 11 月, 米国ソルトレークシティー開催 ) において, ゴードン ベル賞を単独受賞した 3). ゴードン ベル賞ファイナリストには, ピーク性能が 20 ペタフロップス ( 京 の約 2 倍 ) の セコイア ( 米国ローレンス リバモア国立研究所 ) を使い, 同様のダークマターシミュレーションで 14 1 ステップの計算時間 ( 秒 ) , , , , ノード数 図 -4 京 上での GreeM のストロングスケーリング. 横軸が使用したノード数 (= 並列数.8 倍が CPU コア数 ), 縦軸が 1 ステップにかかった計算時間を表す. 破線は理想的なスケーリングが達成された場合の性能である. カラーの違いは評価に用いた粒子数の違いである (2, ,240 3 粒子 ). ペタフロップス ( 実行効率 69%) を達成した米国のグループがあった 4). ところが我々のグループのアプリケーションが実際の計算速度で上回り, ピーク性能が約半分の 京 を用いたにもかかわらず, 1 粒子あたり 2.4 倍の速さでシミュレーションをすることが可能であった ( 米国と同じ計算機を用いた場合は 5 倍近く速い ). こういった点が評価され, ゴードン ベル賞の受賞につながったと考えられる. 我々のアプリケーションが圧倒的に速いのはいくつか理由があるが, 通信部分のアルゴリズムを工夫し並列性能を上げたのに加え, 優れた重力計算アルゴリズムを用いているのが効いている. 米国のグループのアプリケーションでもツリー法を用いている. 我々のものは, ある粒子が受ける力はあわせて数千個程度の質点 ( 粒子と重心 ) からの力として計算すれば, 必要な計算精度としては十分であった. ところが, 米国のグループは, 我々より 5 ~ 6 倍多い数を必要としていた. つまり我々のアプリケーションは, 遠方の粒子群のまとめ方が優れており, そもそも必要な演算量が 1/6 ~ 1/5 で済んでいたのである. 一方, 力の計算の部分は浮動小数点演算の密度が高いため, 実効効率は計算速度とは逆に米国のグループの方が高くなっている. だがこれは無駄な計算をしているとも言える. 単に大規模なスーパーコンピュータ上で実行効率を高めるのではなく, より優れたアルゴリズムを実装し, 実アプリケーション 8 京 の威力で宇宙の正体に迫る ダークマターの超大規模シミュレーション 833

6 特集 続 スーパーコンピュータ 京 の利用 図 -5 京 による世界最大規模のダークマターシミュレーション 左パネルは 現在の宇宙の大規模構造を 5,500 億粒子 を用いてシミュレーションしたものであり およそ 50 億光年にわたる空間の銀河分布を再現している 白枠内は拡大図で 15 あり 右下のパネルは銀河団サイズ 約 10 太陽質量 のハローである 右パネルは 宇宙初期に形成するダークマター構 造を 690 億粒子を用いてシミュレーションしたものである 宇宙誕生からおよそ 1 億年後の宇宙であり 右下のパネルは約 太陽質量のハローである としての性能が優れていたことが評価され 我々の それは対象とするスケールの違いである ハローは 劣勢を覆した勝利につながったと考えられる 密度揺らぎが成長してできるが ダークマター素粒 子が固有に持っている自由運動によって あるスケ 京 による世界最大規模のシミュ レーション 京 コンピュータの非常に強力な演算能力を最 のハローは形成しない ダークマター素粒子として 質量が 100GeV 程度のニュートラリーノを考えた -6 場合 対応するスケールは 10 太陽質量程度である 大限活かすことのできるシミュレーションコードは 左パネルのシミュレーションはそれよりはるかに大 完成した そして我々は 2 種類の世界最大規模の きい銀河や大規模構造の形成を対象としている一方 ダークマターシミュレーションを達成した 図 -5 右パネルのシミュレーションはこのスケールを対象 左パネルは現在の宇宙の大規模構造を 5,500 億粒子 としているので できるハローの数が少なく 構造 を用いてシミュレーションしたものである 右パネ が滑らかになる なお両スケールは質量で 20 桁以 ルは宇宙初期に形成するダークマター構造を 690 億 上の差があるので 両方を同時にシミュレーション 粒子を用いてシミュレーションしたものである い することは エクサスケールをもってしても不可能 ずれも それぞれのスケールのシミュレーションで である は世界最大粒子数である 両者は似ているが いくつか大きな違いがあるこ こういった大規模シミュレーションから何が分か るのだろうか? まず 左パネルのような宇宙の大 とが分かる 宇宙初期に形成するダークマター構造 規模構造形成シミュレーションからは 銀河や銀河 は 現在に比べ滑らかな構造を示していることが分 中心の超巨大ブラックホールの空間分布を予測する かる フィラメント構造の中や ハロー周囲に存在 ことができる 超巨大ブラックホールは銀河と共進 する小さいハローの数が少ないことも分かる 化していると考えられているが その形成 進化過 これらの違いは一体どこからきたのであろうか? 834 ール以下の密度揺らぎは減衰し そのスケール以下 程は謎に包まれている シミュレーションから得ら

7 れたハローの形成史をもとに, 現象論的なモデルを用いて, 個々のハローでできる銀河やブラックホールの性質を理論的に予言することができる. 京 以前のシミュレーションでは, 分解能と体積の両方が足りず, 銀河の階層的形成を追いつつ, 最新の観測データと比較できるような広大な領域にかけて, 銀河とブラックホールの分布両方を予測することができなかった. 今回 京 の威力によりそれがはじめて可能となり, 共進化過程の解明への手がかりなどが得られるようになった. 右パネルのシミュレーションからは, 銀河系が存在するハローの微細構造に関する手がかりが得られる. ここ十数年の研究により, ハローは中心が高密度であり, さらにその中には図 -5 でも見られるような無数のより小さいハローが存在することが分かってきた. また最近の研究は, 小さいハローほど内部質量の中心集中度が大きく, 大きいものと比べて多量に存在することを示唆している. したがってこの高密度な最小のハローが, 地球近傍に存在し得るかどうかが, ニュートラリーノの探査とその詳細な性質を解明するために重要である. なぜなら, ダークマター対消滅起源のガンマ線検出のために, 地球からできるだけ近くて, ダークマターが集中して存在すると考えられる領域に観測衛星を向けたいからである. だが 京 以前は最小のハローのみを正しくシミュレーションできただけであり, 最小のハローがより大きいハローの中でどう分布するかが分からなかった. 今回のシミュレーションでは最小ハローが合体してできるハローを追うことができ, 銀河サイズまで成長したとき, そのなかでの最小ハローの分布をある程度推定することが可能となった. 図 -5 左で示したシミュレーションは, 一辺 50 億光年の非常に広大な領域をシミュレーションしているとはいえ, 宇宙の一部分である. 宇宙のサイズにせまる領域をシミュレーションできるようになることで, 我々の宇宙の構造形成モデルをより精密に決められるようになるだろう. また銀河系のハローのより微細な構造を明らかにすることができ, ダークマター探査のための, より具体的なモデルを提案することができるようになるだろう. またダークマターシミュレーションの規模を上げるだけではなく, バリオンを同時にシミュレーションすることも重要である. 今後は, 今回開発したアプリケーションを, 圧縮性流体のための代表的な計算法である SPH 法 5) などのより複雑な物理を扱うアルゴリズムと組み合わせ, 星形成, 銀河形成など, 宇宙の構造形成過程に関する科学的成果の創出に向けていきたいと考えている. 参考文献 1) 宇宙論 II : 宇宙の進化, シリーズ現代の天文学, 日本評論社 (2007). 2) Barnes, J. and Hut, P. : A Hierarchical O (N log N) forcecalculation Algorithm, Nature, Vol.324, pp (Dec. 1986). 3) Ishiyama, T., Nitadori, K. and Makino, J. : 4.45 pflops Astrophysical N-body Simulation on K Computer : The Gravitational Trillion-body Problem. In Proceedings of the International Conference on High Performance Computing, Networking, Storageand Analysis, SC'12, pp.5:1-5:10, Los Alamitos, CA, USA, IEEE Computer Society Press (2012). 4) Habib, S., Morozov, V., Finkel, H., Pope, A., Heitmann, K., Kumaran, K., Peterka, T., Insley, J., Daniel, D., Fasel, P., Frontiere, N. and Lukíc, Z. : The Universe at Extreme Scale : Multi-petaflop Sky Simulation on the bg/q. In Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis, SC'12, pp.4:1-4:11, Los Alamitos, CA, USA, IEEE Computer Society Press (2012). 5) シミュレーション天文学, シリーズ現代の天文学, 日本評論社 (2007). (2014 年 3 月 1 日受付 ) 8 京 の威力で宇宙の正体に迫る ダークマターの超大規模シミュレーション エクサ時代の構造形成シミュレーション きたるエクサスケールの時代では, 我々は 京 の 100 倍の演算能力を手に入れることができるだろう. 粒子数が 100 倍になれば, さらに色々面白いサイエンスを展開できるようになる. たとえば, 石山智明 ishiyama@ccs.tsukuba.ac.jp 1982 年生.2005 年東京大学教養学部広域科学科広域システム科学系卒.2007 年同大学院修士課程修了.2010 年同大学院博士課程修了. 学術博士. 日本学術振興会特別研究員を経て,2011 年より筑波大学計算科学研究センター研究員. 宇宙の構造形成,HPC の研究に従事.2012 年 ACM Gordon Bell 賞受賞. 天文学会. 835