先端的固体地球シミュレーションコードの開発
Development of Advanced Simulation Tools for Solid Earth Sciences
古市 幹人 (海洋研究開発機構 地球内部ダイナミクス領域)
Mikito Furuichi
このプロジェクトの目的は、数々の固体地球科学の重要な未解決問題に対して、他に類をみない高度なシミュレーション研究を行うために、大規模並列計算機向けの先端的なシミュレーション手法、及びコードを開発することである。これまでに我々は、新しい球面格子「インヤン格子」、マントル対流の新解法「ACuTE法」、粘弾性流体の新しい解法である 「CIP-CSLR-CS法」さらに、局所的粘性差にロバストなストークス流の解法の開発「SSDDソルバー」の開発に成功した。本年度も引き続きこれらのオリジナルの計算手法の発展に取り組む。また、開発した手法を地球科学の具体的課題に応用し、各現象のシミュレーションによる再現に挑戦する。
本年度の主な目標として、前年度までに開発した数値惑星コードの具体的な適応問題として、自由境界・自己重力場での3次元レイリーテイラー不安定性による金属層の落ち込みを再現する。そして、地球型惑星史全体にとって重要な問題である地球核の形成過程の時間スケールに対して新たな制約条件を与えることを試みる。
本年度は、昨年度に開発した、地球の自転速度の時間的変動の効果を取り入れたダイナモモデルを用いてプロダクトランを行う。地球の自転速度には数十年から数万年以上の周期に及ぶ変動が見られる。気候変動による大陸氷床の消長と水の移動が地球の自転速度変動をもたらし、それが核の対流の変動と地球磁場変動を引き起こすという考えが提唱されている。しかしながら、これまでに自転速度変動の効果を考慮したダイナモモデルは無く、この効果を取り入れた世界初の成果となる事を目指す。
マントル対流は、多成分系・固液2相系のマントル対流の3次元シミュレーション向けのプログラムの開発と機能向上、及び中規模モデルによるテスト計算を中心に行う。
当初計画
プレート境界形状・三次元媒質分布を考慮した3次元の地震サイクルシミュレーションコードの開発・ES向けのチューニングを実施する。
修正計画
東北地方太平洋沖地震の発生メカニズム解明に向けて、地震に伴うプレート境界でのすべり・すべり分布を観測データから抽出し、将来的な再現シミュレーションの基礎データを作成する。
年度前半は、レイヤー状にマントル、金属、未分化層を分布させたシミュレーションを行い、自己重力場における3次元レイリーテイラー不安定性による金属層の落ち込みを再現した。また後半では隕石衝突を動的なモデルとして組み込んだ計算を行い、粘性等のパラメータを調整させる事でプロダクトランを行った。また年度を通じて開発したソルバーの効率的な実装や、大規模計算に向いた移流手法の開発等をES2上で行った。
観測されている地磁気変動では、変動周期によって変動幅に違いがある(数万年周期で数十%、約60年周期で約1%)。これを調べるため、自転速度変動周期が磁場散逸時間(地球では2万年に相当)とその1/100(同200年)の場合について計算を行った。前年度よりもさらに長時間(磁場散逸時間の20倍以上、地球では40万年以上に相当)の積分を行い、地磁気変動を調べた。また、自転速度変動の振幅を変化させた場合と、エクマン数(=粘性力/コリオリ力)を変化させた場合に自転速度変動の影響がどのようになるかも調べた。
ES2 の対話型ノードを用いて、ACuTE 法プログラムの動作確認を実施した。
東北地方太平洋沖地震の震源域を対象として、GeoFEMにより海底地形・媒質不均質などを考慮したすべりによる地表変形の応答関数を計算した。
海外からのリモートでの利用であるため、チューニング作業が主であった。
マグマオーシャンからの冷却によって地球内部にレイヤー状の構造が出来ると仮定したモデルにおいて、金属層の中心への3次元レイリーテイラー不安定性による落ち込みをシミュレーションにより再現した。そして、その落ち方が初期地球の内部物性によって大きく異なることを示した。また隕石衝突モデルを導入する事によって、地球が動的に1-10Maといった時間スケールで成長する事を再現した。このような3次元コア形成シミュレーションは世界初である。そして、コア形成過程の隕石衝突による溶融・分化プロセスを逃れた領域が存在する可能性をシミュレーション結果の解析により発見し、それらの定量的な評価を行った。この事は、D”層にその存在が示唆されている始源物質（primordial material）の形成過程についての新しい仮説となりえる研究成果である。
数値手法の面では“移動数値境界流れのシミュレーション”において開発しているVOFを導入する事により、更なる大規模計算に必須となる移流手法の改善を行った。またES2向けにソルバーの効率的な実装方法を開発した。詳細は項目8のチューニングによる成果に記す。
地球自転速度変動の周期が磁場散逸時間(2万年)の場合、コア内の磁気エネルギーとダイポール成分が約30％変動し、その1/100(200年)の周期の場合、ダイポール成分が約1%変動するという結果を得た。この変動幅は観測されている地磁気変動(数万年周期で数十%、約60年周期で1%程度変動)とほぼ傾向の合う値であり、周期によって自転速度変動の効果が異なる事が分かった。また自転速度変動振幅(上のケースでは2%)を1/10にすると周期的な磁場変動傾向が弱くなり、さらにエクマン数を10倍程度上げた計算ではより大きな振幅(6%以上)を与えても変動が起きないことが確認された。つまり自転速度変動の影響は変動振幅のみならずエクマン数にも依存する事が分かった。
多成分系・固液2相系のマントル対流3次元シミュレーションプログラムのうち、(i) 固体の流動を計算する機能の ES2 向け改良、及び(ii) 物質の移動を計算する機能の構築、の2つを集中的に行った。(i) では、ベースとなる ACuTE 法を用いた3次元箱型プログラムのうち、ES2での動作が不完全であった箇所 (非弾性流体近似に基づいた解法、及び F90 非準拠の機能) を修正し、ES2 での動作を確認した。また (ii) では、 (a) FCT 法による移流方程式ソルバの導入、 (b)Darcy 則に基づく固相と液相の相対運動の計算、及び (c) 2相系流体の運動と熱対流運動のカップリング、を行う機能を追加した。また中規模モデルによるテスト計算により、その動作を確認した。
東北地方太平洋沖地震の震源域を対象として、海底地形・媒質不均質などを考慮したすべりによる地表変形の応答関数を計算し、それに基づいて、地震時すべり分布・postseismicなすべり分布を求めた。他の解析結果と比較すると地震時に陸側に近い場所で大きなすべりが発生している結果となった。
表面張力モデルを従来よりも高密度側に張力をシフトして計算する事により、計算を安定化させる事に成功した。
(年度当初の研究計画を全て達成した場合を100%として数値で示してください。複数の目標があった場合は、それぞれについて達成度を数値で示してください。)
数値惑星シミュレーション : 90％
地球ダイナモシミュレーション : 80％
マントル対流シミュレーション : 70%
地震発生サイクルシミュレーション : 計画の見直しを行ったため、当初の計画に対しては10%、見直し後の計画に対しては50%
移動数値境界流れのシミュレーション : 50％
年度使用量計画15000ノード時間に対して、2012/01/04時点での使用は全体で約7400ノード時間である。3月末までに計画通りの資源を使用する予定である。
年度使用量計画15000ノード時間に対して、2012/01/04時点での使用は8770ノード時間である。3月末までに計画通りの資源を使用する予定である。
現時点ではES2の対話型ノードの使用に留まっている。
本年は、計画の変更のため、ESの大規模利用はすべり応答計算にのみとなった。ただし、将来的な東北地震の再現に向け、M7〜M9を対象としたマルチスケールサイクルシミュレーションの実現のため、サイクルシミュレーションの計算アルゴリズムを変更し、ESでのチューニングを実施している。
(ベクトル化、並列化チューニング等、計算機資源を有効利用するために行ったこととその効果を記載してください。)
年度初頭にベクトル化、並列化チューニングは終わっており、64ノード申請も通っている。本年度はさらにマルチグリッド法の実装について、グリッドレベルに合わせて実行するMPIコミュニケータを変更する調整を行い、論文にて発表している。さらに、“移動数値境界流れのシミュレーション”と共同でred-black SOR法のループ計算でメモリアクセスを連続にさせるreordering法の開発を行い、現在論文投稿中である。
地球ダイナモシミュレーションコードは、初代地球シミュレータ向けにほぼ限界までチューニングし、2004年ゴードンベル賞を受賞したコードである。それから基本的には何も変更していない。現地球シミュレータでのベクトル化率は99.34％、64ノードと16ノードで計測した並列化率はsuper linear（100％以上）である。
ACuTE 法の根幹である、多重格子法による流れ場解法に含まれるループのベクトル化を実施した。
マルチスケール地震サイクルシミュレーションの実現に向けH-matrix法を導入した。この手法によると演算数・必要メモリをオリジナルのO(N^2)からO(N)〜O(NlogN)程度に軽減できる。しかし、アルゴリズムの特性上、最内ループが短くなるためESでの性能はまだ低い。現状、ESに最適化されたBLASルーチンを使用することにより、約2倍高速化したところである。
並列性に問題のあった接触角の部分を並列化した。具体的には、接触角の条件を満たすように界面をGhost Cellに外挿して、9色のマルチカラーを使う事により並列化した。
物理過程モデルの導入と、パラメータセッティングに手間取ったため、そのベンチマーク等を含めて大規模計算での使用は短時間なものにとどまっているが、それらの問題はおおむね解決したため、今後は計画に沿った運用が期待できる。
年度前半は計算結果の解析等に時間を取られてジョブ投入ペースがやや鈍ったが、現在は様々なエクマン数と変動振幅に対する自転速度変動の効果を調べ、スケーリング則を確立するため大規模パラメータサーベイを行う段階に入っており、2012年1月～3月の3ヶ月間で計画通りの資源を使用する予定である。
プログラムに新たに導入した物理素過程の解法ルーチンの動作を検証する作業を集中的に実施したこと、及び ES2 での動作が不完全であった原因の追求に時間を要したことが原因である。しかしこの結果、問題点の洗い出しがほぼ完了したため、今後の作業は円滑に進むと期待できる。
「愛媛大学『研究室からこんにちは』」(アトラス出版) に紹介記事。