先端的固体地球シミュレーションコードの開発
Development of Advanced Simulation Tools for Solid Earth Sciences
古市 幹人 (海洋研究開発機構 地球内部ダイナミクス領域)
Mikito Furuichi
本プロジェクトの目的は、数々の固体地球科学の重要な未解決問題に対して、独自の大規模並列計算機向けの先端的なシミュレーションコードを開発することで、他に類を見ない計算機性能のアドバンテージを得た高度なシミュレーション研究を行うことである。これまでに我々は新しい球面格子「インヤン格子」、マントル対流の新解法「ACuTE法」、粘弾性流体の新しい解法である 「CIP-CSLR-CS法」さらに、局所的粘性差にロバストなストークス流の解法「SSDDソルバー」といった手法を開発した。これらによって、我々のシミュレーションコードは、従来の研究では扱う事の出来なかった解像度、粘性コントラスト等をES2上で扱う事が出来る。本年度も引き続きこれらのオリジナルの計算手法の発展に共通して取り組んだうえで、開発したコードによる地球科学の具体的課題への応用に挑戦する。
3次元コア形成シミュレーションに状態変数依存の物理モデルを導入し、そのプロダクトラン行う。多成分系・固液2相系を扱う混相流の研究では、マントル対流の3次元シミュレーションコードの開発と機能向上、及び ES2 向けの最適化作業を中心に行う。さらに部分溶融対の発生・移動を伴うマントル対流問題を対象にしたテスト計算を実施する。また、移動数値境界流れのシミュレーションは、計算時間の８割以上を消費するMGスムーサーを改良し、テストする。
昨年度開発した自転速度変動を考慮したモデルにおいて、自転速度変動が地磁気変動を引き起こす事が分かった。変動周期は約2万年（気候変動と関係がある周期であるミランコビッチ・サイクルの1つに非常に近い値である）としたところ数十％の地磁気変動が生じた。これは実際に古地磁気学でも明らかになっている地磁気変動量である。そこで本年度は地磁気変動と自転速度変動がどのように結びついているのかを詳細に解析すると共に、自転速度変動の周期や振幅、エクマン数(=粘性力/コリオリ力)を変化させた場合についても調べていく。
状態変数に依存しない様々なコア形成シミュレーションのセットアップに対して様々なベンチマークテストを、最大1024^3のグリッドサイズで64ノードを用いて行った。一方で温度場を解いて物性を状態変数に依存させた計算もおこなった。その際、一様な物性の場合と比較して問題の収束率が格段に悪くなるため、128^3のグリッドサイズで8ノード計算を実行している。また今年度は粒子―流体混相流コードを新規に開発し、現在8-16ノードを使用して、100～500万粒子の計算を64^3～128^8グリッドの流体にカップルさせた計算を行っている。年度末までには1000万粒子を扱う計算を行う予定である。また、これまでに開発した「ACuTE法」によるマントル対流シミュレーションを改良し、1) 太陽系外惑星スーパーアースのマントル対流のシミュレーション、2) プレートの剛体的運動や局所的破壊、沈み込みを取り扱うためのプログラムの高度化に新規に着手した。(1)については1024*256の解像度で1パラメータあたり2ノードで約150時間以上の計算を15パラメータ程度行っている。 2) については、流体粒子が過去に受けた変形の履歴を記憶する「ダメージパラメータ」を導入するためのモデル開発を行い、256*256*128程度の解像度で8ノードでテスト計算を行っている。解像度を上げる必要があるため、さらなる高解像度計算を行う予定である。
自転速度変動を考慮した地球ダイナモシミュレーションについて、地磁気が変動するメカニズムを調べた。その過程でより詳細なデータや当初着目していなかった量のデータを取る必要が生じたため追加の計算を実施した。その結果、地磁気以外にも変動する重要な物理量の存在、地磁気変動が自転速度変動によりどのように生じるかのメカニズム、地磁気変動とそれ以外の変動量との関係などについて明らかにする事が出来た。自転速度変動を考慮したダイナモモデル自体が世界で最初のものであるため、これらの発見は画期的な成果となると考えている。さらに自転速度変動の周期を変化させた場合やエクマン数を変えた場合についても計算を行った。
本年度は開発したコア形成シミュレーションコードの解像度依存性などのベンチマークテストを行うと同時に、様々なインパクトサイズ、プロトプラネットの粘性強度といったパラメータに対してシミュレーションを行い、集積プロセスにおいて火星サイズのプロトプラネットがどの程度存在できたのかを定量的に評価することを試みた。その結果、N体計算から示唆される幅広いパラメータ範囲において、我々が主張したい月形成に資するジャイアントインパクト以前に隕石衝突による直接的な溶融・分化プロセスを逃れた相当量のプロトプラネットの残骸が存在しうることを見出した。さらに、エネルギー方程式を導入することで、コア物質の地球中心への落ち込み時の、温度状態の変化を取り扱うことを可能にした。さらにMGスムーザーをSOR型に取り換え、高実行性能での使用が可能なように改良した。
また、現在の地球内部の物質分化においてキーとなる深部マグマだまりや、ジャイアントインパクトに伴うマグマオーシャンでの溶融・固化プロセスのダイナミクスを理解するために、これまで開発してきた浸透流コードに加えて、今年度新たに粒子―流体混相流コード開発を行った。具体的には粒子部分を、多粒子を高速に取り扱える離散化要素法DEMによって解き、流体部分に高粘性で有効なストークス流を用いてカップルさせた。流体―粒子の相互作用には現象論に基づく平均化モデルを使用した。このような混相流コードの研究は工学分野において長年行われてきたが、地球科学において関心のある高粘性流体に対して行われた例はなく、従来の手法では高粘性領域を解くことが困難であるため、新しく数値手法開発を行った。高粘性を扱うための工夫として[1]DEMの加速度項の除去[2]粒子に作用する流体圧力への半陰解法の適応[3]バルクの流体粘性モデルと高粒子密度領域における有効粘性モデルとのハイブリッド解法を実現することで、[1]高粘性中でのDEM粒子の求解[2]ソルバーの安定性強化[3]粒子―流体間の同期コストの削減を行うことに成功した。
さらに本プロジェクトおいてこれまでに開発した「ACuTE法」によるマントル対流シミュレーションコードを改良し、1) スーパー地球におけるマントル対流の数値シミュレーション、2) プレートの剛体的運動や局所的破壊、沈み込みを取り扱うためのシミュレーションモデルの開発にも取り組んだ。
1) については、太陽系外惑星で見つかっている最大で質量が地球の10倍程度もある「スーパー地球」のマントル対流がどのようになるかを調べた。マントル対流はその惑星の表層環境を強く規定することから、系外惑星のマントル対流を調べることは「その系外惑星がhabitable (生命体が居住可能) な環境であるか否か」という究極的な問題につながるものである。ここではACuTE法の得意とする高レイリー数及び粘性の強い温度依存性を持つ流体の熱対流に、スーパー地球のマントル物質に期待される強い圧縮性を考慮に入れるように改良した手法を用いて2次元モデルによるシミュレーションを行った。これまでレイリー数や粘性の温度依存性を様々に変えた15ケースほどの計算を行っている。レイリー数は最大で10の10乗、質量はスーパー地球でほぼ最大である10倍地球質量とした。これらは扱いが困難であり他の研究例が未だ存在しない領域である。結果の詳細は現在解析中であるが、スーパー地球のマントル対流では地球のそれと比べて熱輸送効率がかなり低い（地球の1/3以下程度）事が分かった。これは惑星内部が極めて冷えにくい事を意味し、スーパー地球の熱的進化の過程に重大な影響を及ぼすと考えられる。また、マントル内に対流がほとんど生じていない「成層圏」のような領域が生じているらしい事も分かった。これは地球のマントル対流とはかなり異なった特徴である。これらの成果は本年度内または来年度早々に論文にまとめる予定である。
2) については、プレートの剛体的運動、局所的な破壊による沈み込み帯の発生などを適切に扱うべく、流体粒子が過去に受けた変形の履歴を記憶する「ダメージパラメータ」をACuTE法に導入した。このモデルの究極の目的は、沈み込み帯からマントル内に取り込まれた水の影響 (マントル物質の粘性率低下、融解、マグマ発生など) を考慮したマントル対流のシミュレーションを行うことであり、「ダメージパラメータ」の導入はその最初の段階としての意味を持つ。さらにこのプログラムを用いて予備的な3次元計算を行った。本項執筆時点では熱拡散時間の0.1% (地球では約2億年に相当) 程度までの計算であるが、表層の温度が低く固いプレート面にダメージが蓄積し、また上昇流付近に強いダメージが集中するなど、テスト計算が順調に進んでいる。ただし現在までのところは「ダメージ」の蓄積に伴うプレートの局所的な破壊 (粘性率の低下) は顕著には見られていないが、さらなる長時間計算を実施し、計算が安定に行えるかを確認していく。その後実際に地球で見られるプレート運動やプレート同士の衝突による沈み込み帯形成などの研究に応用していきたいと考えている。
自転速度変動により生じる事が分かった地磁気(地球コア内の磁気エネルギー及び地磁気のダイポール成分)の変動について、変動が生じるメカニズムを調べるため詳細なデータを取得し結果を解析した。その結果明らかになった事について要点を述べる。(1)地磁気変動に伴いコア―マントル境界の熱流量も変動する事が分かった。変動振幅は10%以上とかなり大きい。(2)自転速度の変動はまずコア内の流れに影響を与える。その結果対流の運動エネルギーも変動する。 (3)自転速度変動と地磁気変動の位相は一致しない(すなわち自転速度が最大の時に地磁気強度も最大といったような単純な関係ではない)。位相差は約π/2である(地磁気強度の最大/最小は平均的な自転速度の時間である)。自転速度の時間微分との関係をとると、これは運動方程式中に新たに現れる重要な量でもあるが、よりはっきりした関係が見られる。位相差は約πである（地磁気強度が最大時に自転速度微分が最小になる）。自転速度と地磁気の変動位相が一致しないという事実は、古地磁気学における地磁気変動の解析にも示唆を与えられるかもしれない。(4)ダイナモ活動やジュール散逸といった量も自転速度変動に伴い変動する事が分かった。これらは磁気エネルギーを変動させる要因であり、解析により磁気エネルギー変動との関係を明らかにする事が出来た。 (5)地磁気ダイポール成分が弱くなった時は、高次モード成分の相対的強度が強くなる事が分かった。これは実際に地磁気エクスカーションでも見られている特徴である。
以上に加えて、より短周期の変動やエクマン数を変えた場合の計算結果も得られており現在結果の解析中である。
(年度当初の研究計画を全て達成した場合を100% / 複数の目標があった場合は，それぞれについて達成度を数値で記載)
数値惑星シミュレーション : 80％
地球ダイナモシミュレーション : 80%
当初計画している年度使用枠15,000時間内で約10,000時間を使用している。年度末に向けて大規模ジョブが走るため、計画に沿った利用が行えていると考える。
年度使用予定枠25000時間の内、本項執筆時点で約12500時間を使用している。1月～3月期にも大規模ジョブを計画しているため、予定通り全時間を使用する見込みである。
(ベクトル化、 並列化チューニング等、 計算機資源を有効利用するために行ったこととその効果を記載)
コア形成シミュレーションコードでは、連続メモリアクセスを可能にするSOR法の実装、マルチグリッドでの荒い解像度におけるMPIプロセス数の最適化やPICやVOFなどの数値移動境界手法の並列化などを実装することでES２でのチューニングを行った。これらのことにより、たとえば8ノード計算において14％の性能が出ている。
一方で、粒子―流体コードの開発においては西浦らが開発したPGS法によるDEMをES２向けにチューニングを行った。その結果ベクトル化率９９．９％でのDEMの実装が行えている。
該当なし．