18:15 〜 19:30
[STT60-P02] 大規模・高詳細な準動的地震サイクルシミュレーションに向けた高並列化の検討
キーワード：地震サイクル, ケイパビリティコンピューティング, 並列計算, 階層行列
Ohtani et al.(2011)は，準動的な地震サイクルの問題に，H-matrices 法と呼ばれる，密行列を階層的な小行列に分割し個々の小行列を効率的に圧縮する手法を適用した．これにより，断層の離散化数Nが105-106の範囲で，M8の地震サイクルを計算した場合の演算数が，従来法のO(N2) からO(N)?O(NlogN) に減少できることとなり，京コンピュータなどの大規模並列計算機を利用したcapacity computingによって，様々なM8クラスの地震発生シナリオを評価できるようになってきた．
しかし，実際のプレート境界では，蓄積される歪(すべり欠損) は本質的に規模が異る地震同士の相互作用によって解消されており，地震サイクルのシミュレーションを現実に近づけるにはこういった様々なスケールの地震の相互作用をモデル化する必要がある．そのためには，対象とする最小規模の地震を解像するだけの細かな空間離散化をモデル化するプレート境界に適用する必要が生じ，例えば，従来の研究で対象としていた地震のマグニチュードを基準とし，そのマグニチュードより2小さな地震まで含めたシミュレーションを実施するとすれば，モデル地震断層の断層長が従来の研究の凡そ1/10になり，プレート境界を100 倍細かなメッシュへ離散化することが要求される．そういった計算は,もはやシリアルCPUでの計算は不可能となり，京コンピュータ等の大規模並列計算機の大部分を利用するような計算(capability computing)を実施する必要が生じる．以上から，現実的な地震サイクルシミュレーションの実現には大規模並列化が必要不可欠である．
我々はこれまで，Ohtani et al.(2011)で扱われたモデルと同程度のモデル断層(N=3x105)にHmatrices 適用し，並列数を百程度とした数値計算を実施してきている．その際，Hmatrices 化前の密行列を基準に行方向に密行列を均等に分割するようなバンド状の部分領域を設定し，各部分領域をMPI プロセスに割り当て，そのバンド幅にオーバーラップする小行列をそのMPI プロセスに割り当てるような一次元分割で並列計算を実施している．この並列化方法は複数のバンド領域を跨ぐ小行列に関するバンド領域間通信が不要となり，コーディングが単純化できる反面，こういった重複小行列に関する演算の一部は，複数バンド領域で同一の計算を冗長に繰り返す必要が生じる．つまり，大並列を仮定した場合，シリアル計算と比較してトータル演算量が著しく増大してしまいストロングスケーリングが成り立ちにくくなるといったデメリットがあり，現行の並列化方法は大並列での実行に適さない．
このことから，今回我々は，行方向のバンド分割をある分割数に抑えることにより重複演算の増加を抑制し，各行バンドに対し，列方向へも領域分割を適用し，分割を2 次元化する方針で並列化を見直した．つまり，従来の一次元分割の各行バンド内の演算を並列処理することによって，全体としてのスピードアップを計ろうとしている。
各行バンド内で，小行列の大きさと独立したサイズ(ブロック数) を基準にしてcyclic な分割を行えば，大きな小行列は複数プロセスに分割され，演算の負荷バランスを解消できる．ただし，これによって従来の一次元分割計算に比べると，小行列内での通信と(最終的な行列-ベクトル積の結果を得るための) 小行列間の通信の両方が各行バンドに追加されることなる．しかし，行分割による演算数増加の抑制と，演算を行バンド内で並列処理できること，との兼ね合いによって，現時点でも，同じ並列数(1024MPI プロセス) の計算(N=1.3x106)に対し，従来法より2 倍以上の高速化を達成している．
今回採用した並列化方法は，ブロック数・各次元方向の並列数等の並列パラメタの指定により演算-通信のバランスが変化し，全体としての計算性能に影響を及ぼすと考えられる．今後は，並列化の更なる効率化と，最適な並列パラメタの探索を実施していくことになる．
謝辞．本研究の一部には理化学研究所の京コンピュータを使用させて頂きました(課題番号hp120278)．地震サイクルコードの並列化・チューニングには富士通(株)のチューニングチームに助力して頂きました．
しかし，実際のプレート境界では，蓄積される歪(すべり欠損) は本質的に規模が異る地震同士の相互作用によって解消されており，地震サイクルのシミュレーションを現実に近づけるにはこういった様々なスケールの地震の相互作用をモデル化する必要がある．そのためには，対象とする最小規模の地震を解像するだけの細かな空間離散化をモデル化するプレート境界に適用する必要が生じ，例えば，従来の研究で対象としていた地震のマグニチュードを基準とし，そのマグニチュードより2小さな地震まで含めたシミュレーションを実施するとすれば，モデル地震断層の断層長が従来の研究の凡そ1/10になり，プレート境界を100 倍細かなメッシュへ離散化することが要求される．そういった計算は,もはやシリアルCPUでの計算は不可能となり，京コンピュータ等の大規模並列計算機の大部分を利用するような計算(capability computing)を実施する必要が生じる．以上から，現実的な地震サイクルシミュレーションの実現には大規模並列化が必要不可欠である．
我々はこれまで，Ohtani et al.(2011)で扱われたモデルと同程度のモデル断層(N=3x105)にHmatrices 適用し，並列数を百程度とした数値計算を実施してきている．その際，Hmatrices 化前の密行列を基準に行方向に密行列を均等に分割するようなバンド状の部分領域を設定し，各部分領域をMPI プロセスに割り当て，そのバンド幅にオーバーラップする小行列をそのMPI プロセスに割り当てるような一次元分割で並列計算を実施している．この並列化方法は複数のバンド領域を跨ぐ小行列に関するバンド領域間通信が不要となり，コーディングが単純化できる反面，こういった重複小行列に関する演算の一部は，複数バンド領域で同一の計算を冗長に繰り返す必要が生じる．つまり，大並列を仮定した場合，シリアル計算と比較してトータル演算量が著しく増大してしまいストロングスケーリングが成り立ちにくくなるといったデメリットがあり，現行の並列化方法は大並列での実行に適さない．
このことから，今回我々は，行方向のバンド分割をある分割数に抑えることにより重複演算の増加を抑制し，各行バンドに対し，列方向へも領域分割を適用し，分割を2 次元化する方針で並列化を見直した．つまり，従来の一次元分割の各行バンド内の演算を並列処理することによって，全体としてのスピードアップを計ろうとしている。
各行バンド内で，小行列の大きさと独立したサイズ(ブロック数) を基準にしてcyclic な分割を行えば，大きな小行列は複数プロセスに分割され，演算の負荷バランスを解消できる．ただし，これによって従来の一次元分割計算に比べると，小行列内での通信と(最終的な行列-ベクトル積の結果を得るための) 小行列間の通信の両方が各行バンドに追加されることなる．しかし，行分割による演算数増加の抑制と，演算を行バンド内で並列処理できること，との兼ね合いによって，現時点でも，同じ並列数(1024MPI プロセス) の計算(N=1.3x106)に対し，従来法より2 倍以上の高速化を達成している．
今回採用した並列化方法は，ブロック数・各次元方向の並列数等の並列パラメタの指定により演算-通信のバランスが変化し，全体としての計算性能に影響を及ぼすと考えられる．今後は，並列化の更なる効率化と，最適な並列パラメタの探索を実施していくことになる．
謝辞．本研究の一部には理化学研究所の京コンピュータを使用させて頂きました(課題番号hp120278)．地震サイクルコードの並列化・チューニングには富士通(株)のチューニングチームに助力して頂きました．