世界最高精度での全球地震波伝播シミュレーションに成功
～理論地震波形を用いた地球内部構造研究が大きく進展～
研究概要地球内部構造決定の高精度化のためには、地震発生により地球全体を伝播する地震波を理論的にシミュレーションすることが必要です。地球はもし完全な球であれば、その上を伝わる地震波は理論的に表すことができますが、実際にはわずかに球からずれているために理論的に表すことはできず、コンピュータによる近似的な解き方により理論地震波形をシミュレーションすることが一般的です。
研究グループでは2003年に当時世界最速のスーパーコンピュータであった「地球シミュレータ」を用いて、初めて完全な球ではない現実的な地球モデルに対して地震波形を周期5秒の精度で全球地震波伝播シミュレーションすることに成功しました。しかしながら、地球上を伝播する地震波の解析には周期1秒の実体波(P波およびS波)が必要なことが知られており、周期1秒の精度での理論地震波形記録シミュレーションは大きな目標でした。
今回、スーパーコンピュータ「京」を用い、理論地震波形記録シミュレーションのための数値解法であるスペクトル要素法のプログラムを最適化することにより、周期約1.2秒の精度という世界最高精度での全球地震波伝播シミュレーションに成功しました。
この成果は、周期1秒の精度で全球地震波伝播シミュレーションを実行することがコンピュータの性能向上により近い将来、実現可能であることを示しており、地震発生メカニズムの解明および地球内部構造の精密な決定の応用への可能性を示した点で、地震や津波などの防災・減災への大きな貢献が期待されます。
本成果は、計算機科学の専門誌「International Journal of High Performance Computing Applications」に2月29日付け（日本時間）でオンライン掲載されました。
背景地球を構成する岩石の弾性的性質により、地震波のP波およびS波では地震波形データを解析するための地震波形は周期1秒の精度が必要であることが知られています。したがって、観測された地震波形を再現する理論地震波形をシミュレーションする上で周期1秒の精度は到達すべき大きな目標でした。近年、大型計算機の発展と共に、コンピュータによる近似的な解法を使用して理論的な地震波形をシミュレーションする手法は大きな進歩を成し遂げてきました。研究グループでは2003年に当時世界最速のスーパーコンピュータであった「地球シミュレータ」を用いて、初めて完全な球ではない現実的な地球モデルに対して周期5秒の精度で全球地震波伝播シミュレーションすることに成功しました(平成15年12月5日既報)。全球地震波伝播シミュレーションでスペクトル要素法によるシミュレーションを実施する際、地球モデルの分割は地球全体を6個の四角錐に分割し、それぞれの四角錐をさらに細かい四角錐に分割してスーパーコンピュータの個々のCPUに割り当てて計算を実行します（図1）。
現実的な地球モデルに対して全球を伝播する地震波形を周期1秒の精度でシミュレーションするためには、このような分割を2003年当時の54億個の格子点に対して、さらに数十倍の細かさで行うことが必要であり、計算規模が大きすぎるために莫大なCPU数とメモリが必要でした。
図1. スペクトル要素法による地球モデルの分割。スペクトル要素法では、まず地球全体を6個の四角錐に分割し、それぞれの四角錐をさらに細かい四角錐に分割しスーパーコンピュータの個々のCPUに割り当てて計算を実行する。2003年に地球シミュレータ上で実施した図の例では、6個の四角錐をさらに26×26＝676個に分割し、全体では676×6= 4056個のCPUを用いた大規模計算を実行した。
成果今回のシミュレーションでは、「京」の82,134ノード(全ノードの99%)を用い、「地球シミュレータ」で理論地震波形記録シミュレーションを実行してきたスペクトル要素法のシミュレーションプログラムを用いて、現実的な地球モデルを6,652億個の格子点に分割することで周期約1.2秒の精度での全球地震波伝播シミュレーションに成功しました。これまでの記録は2008年に米国のグループが達成した約2秒の精度でのシミュレーションでした。
「京」上で大規模なシミュレーションを実行して十分な性能を発揮するためには、スペクトル要素法のプログラムを最適化することが必要でした。今回、プログラム中の計算量の8割を占める部分の最適化を図り、「京」の1ノードである8CPUを用いて計算量を分割した場合の性能向上が7.89倍となりました(図2)。さらに、並列計算での性能を示すストロング・スケーリングでも、同じ規模の計算を36,504ノードで実行した場合と比較して、82,134ノードでの結果は99.54%の性能を示すという、極めて良好な結果が得られました。2003年当時、地球シミュレータの243ノードを用いたシミュレーションの実効性能は5テラフロップスでしたが、82,134ノードを用いた計算では、実効性能が1.24ペタフロップス(ピーク性能比11.84%)という高い性能を実現しました。
このシミュレーションでは、2011年東北地方太平洋沖地震(マグニチュード9.0)で生じた地震波を日本列島付近の地震観測点で観測された地震波形と比較することを試みました(図3)。長野県松代の地震観測点における地震波形と理論波形を比較した結果は良い一致を示しており、地震の初期破壊の様子を周期1.2秒の精度でシミュレーションした理論地震波形が良く再現していることを示しています。
図2. スペクトル要素法のプログラムの最適化の結果。青で示した結果はオリジナルのコードに対し「京」のコンパイラによる自動並列化のみを適用した場合の結果であり、当該部分の処理がほとんど並列化されていない。赤で示した結果はプログラムを最適化して効率的な並列化を行ったため、当該プログラム部分の並列化が可能となった。
図3. 長野県松代の地震観測点で観測された2011年東北地方太平洋沖地震の発生により生じた地震波と理論地震波形記録シミュレーションの比較。黒線が観測波形、赤線が理論地震波形を示す。波形は地動の上下動速度成分(単位はm/sec)を示し、地震発生時から約3分間の記録を示している。地震の破壊過程を反映した観測波形の山と谷の様子が周期1.2秒の精度で再現されている。
今後の展望今後、周期1秒の精度でのシミュレーションを実現するためには少なくとも、110,000ノードでの実行が必要と推定されますが、「フラッグシップ2020プロジェクト（ポスト「京」の開発）」で進められている次世代のスーパーコンピュータでは十分実現可能な計算資源と考えられます。今回の結果に基づいて、今後も周期1秒の精度での理論地震波形記録シミュレーションの実現に向けてシミュレーションの高度化を進め、地震発生メカニズムの解明および地球内部構造の精密な決定を通して、地震や津波などの防災・減災へ貢献していきたいと考えています。
)