令和３年度高性能汎用計算機高度利用事業「富岳」成果創出加速プログラム「大規模数値シミュレーションによる地震発生から地震動・地盤増
幅評価までの統合的予測システムの構築とその社会実装」
成果報告書

令和４年５月３０日
国立研究開発法人海洋研究開発機構
堀高峰

本報告書は、文部科学省の高性能汎用計算機高度利用事業費補助金事業として、国立研究開発法人海洋研
究開発機構が実施している「富岳」成果創出加速プロ
グラム「大規模数値シミュレーションによる地震発生
から地震動・地盤増幅評価までの統合的予測システムの構築とその社会実装」の令和３年度の成果を取りまとめたものです。

目次
１．補助事業の目的 .................................................................... - 1 ２．令和３年度（報告年度）の実施内容 .................................................. - 1 ２－１．当該年度（令和３年度）の事業実施計画 ........................................................................... - 1 （１）国の被害想定のための地震・津波災害予測でのアプリケーション活用（２）地震に関する災害被害予測のための大規模アプリケーションの開発（３）開発した計算アプリケーション群を実務で活用するための環境整備（４）プロジェクトの総合的推進
２－２．実施内容（成果） ............................................................................................................... - 2 （１）国の被害想定のための地震・津波災害予測でのアプリケーション活用（２）地震に関する災害被害予測のための大規模アプリケーションの開発（３）開発した計算アプリケーション群を実務で活用するための環境整備（４）プロジェクトの総合的推進
２－３．活動（研究会の活動等） .................................................................................................- 115 ２－４．実施体制- 115 -

補助事業の名称「富岳」成果創出加速プログラム「大規模数値シミュレーションによる地震発生から地震動・地盤増幅評価までの統合的予測システムの構築とその社会実装」

１．補助事業の目的
本研究では、国難とされる首都直下地震や南海トラフ地震に備えるために必要となる、地震災害の定
量的評価手法を、国が実施する地震動・地盤増幅評価に実装することを目指す。そのために、スーパーコンピュータ「富岳」での超大規模計算を実施するための統合的予測システムを構築する。それとともに、地震発生予測に工学的な品質保証の概念を導入して政策判断等に耐える基礎を築く。そのうえで、
従来、重要構造物の設計に際し、同じ想定地震に対して個別に地震動・地盤増幅の評価を行ってきた建築・土木系の企業が、国の被害想定のために実装・使用されるものと同等の計算ツールやデータを実務で活用するための仕組みを構築する。

２．令和３年度（報告年度）の実施内容
２－１．当該年度（令和３年度）の事業実施計画（１）国の被害想定のための地震・津波災害予測でのアプリケーション活用
重点課題から開発してきた大規模有限要素コードを用いた長周期地震動の計算結果が、国の被害想定に向けた委員会でこれまで検討対象となってきた。本年度は、国の被害想定に向けた委員会からの長周期地震動計算結果の公表を目指し、同様な計算での精度検証作業やそのために必要な計算アプリケーシ
ョンの修正を行う。
「富岳」を用いた計算が予定される非線形地盤増幅評価に関しては、評価に対して必要となるより高分解能な地震動計算を、大規模有限要素コードを用いて試算するとともにポスト処理を行う。さらに、弾性波動場ならびに弾性地殻変動応答の参照解の計算を実施する。
（２）地震に関する災害被害予測のための大規模アプリケーションの開発「富岳」の性能を引き出すように、計算科学・計算機科学の最先端技術を駆使することで、地震に関する災害被害予測のためのアプリケーション群の基幹数値解析コンポーネントである、地震・地震動・
地殻変動計算等を行う大規模数値解析用コード群に対し、「富岳」実機上でより良好な性能を発揮するよ
うに昨年度に引き続き研究開発を行う。本年度は、「富岳」の計算機構の特性を考慮した上で、「富岳」の性能を引き出すように、solver の前処理等の改良を行い、アプリケーション群の性能向上を図る。また、超大規模解析モデルを構築し、これを用いた解析を行うことで、解析手法の有効性を示す。
（３）開発した計算アプリケーション群を実務で活用するための環境整備
国の被害想定算出に使われるものと同等の長周期地震動計算・非線形地盤増幅計算を、建築・土木系の各企業が独自の問題設定で実施するためのサポート体制を、海洋研究開発機構が中心となって整備する。昨年度、計算アプリケーション群を使用するための環境構築をすでにある程度まで完了し、プロジ
ェクトに連携機関として入っていただいている企業の方々にベンチマーク問題を対象に試用していただ
き、フィードバックをいただいている。本年度は、フィードバックをもとに、構築した計算環境のさら-1-なる改良を行う。またより現実的な地下構造を対象とした実務に近い問題への適用を企業の方々と連携して進めることを通じ、さらなる機能の充実を図る。
（４）プロジェクトの総合的推進
プロジェクトの総合的推進にあたっては、参画各機関の協力・連携を密にし、円滑にプロジェクトを運営していくため、進捗管理等のための会合等を開催するとともに、外部有識者を招聘した諮問委員会
から、本プロジェクトの進捗状況の把握・評価・改善提言・指導等を受ける。また、プロジェクトで得られた成果については、広報や研究活動を通じて積極的に公表するとともに、「富岳」成果創出加速プロ
グラム事業全体のアウトリーチにも協力する。さらに、重点課題で開発したアプリケーション群や各関係機関で得られた貴重なデータならびに、本年度の業務で得られる成果等を保管・管理し、研究を推進するための計算機環境を維持・整備する。
２－２．実施内容（成果）（１）国の被害想定のための地震・津波災害予測でのアプリケーション活用
本プロジェクト参加者らは、東大地震研のグループを中心に、先行課題であるポスト「京」重点課題から継続して、「富岳」の計算性能を最大限活用できるような大規模有限要素計算アプリケーションの開発に取り組んできた。重点課題において開発された地震動解析プログラム、非線形地盤増幅計算プログ
ラム、地殻変動計算プログラムは、国の被害想定のためのハザード評価における計算手法高度化のニーズと合致しており、国の被害想定への適用を目指した検討にすでに組み込まれている。これらに基づき、「富岳」でしかできない規模の地震動・地盤増幅・地殻変動計算による被害予測システムを構築し、国の被害想定のためのハザード評価の過程に実装することが本課題における大きな目標となる。本項目では、まず①において、国の被害想定のためのハザード評価、具体的には長周期地震動評価に対する本プ
ロジェクトの協力について紹介する。②において、①のための長周期地震動評価に使用する上で必要な精度検証計算を行った内容を述べる。③において、津波評価や地震シナリオのベースとなる地殻変動計
算について、数値誤差の定量化を行うための、「富岳」での大規模計算に基づいた検討について述べる。
①内閣府防災の「相模トラフ沿いの巨大地震に伴う長周期地震動による被害予測の高度化」への協力
将来の発生が懸念される海溝型巨大地震による長周期地震動の推定は、都市部の地震被害対策を進めていくうえで重要である。そのために、適切に設定された震源断層モデルに対し複雑な地盤ならびに地殻構造における地震波動場を数値計算することが必要である。国の被害想定では従来、差分法による地震動解析手法が使われてきた。一方、差分法は構造格子に基づいて計算を行うため、速度コントラストのある層境界面の形状を忠実に再現できず、結果的に計算精度が悪くなる場合があることが指摘されている[1]。このような場合、有限要素法のような非構造メッシュに基づいた計算手法を用いることが望ましい。ここで、有限要素法導入の障壁となっていたのは、計算コストの大きさだった。本プロジェクトの東京大学地震研究所のグループが開発してきた地震動の大規模有限要素計算用プログラム (例えば[2][3][4])は、四面体二次要素による詳細な形状表現、陰的時間積分による安定な計算、高性能前処理付き
反復解法による高速かつスケーラブルな一次方程式求解により、複雑地盤における大規模地震動計算を精度よく、安定に、かつ高速に行うことを可能としている。その高性能な solver はこれまで国内外のフ-2-

ラグシップレベルのスーパーコンピュータにおいて良好な性能を出しており、高性能計算の分野で高い評価を得ている。また、同グループは、成層構造地盤に最適化した全自動大規模有限要素メッシュ生成
ツールも開発しており、計算モデルの生成に関わる計算コストの問題も解決している。そのため、国が行うような被害想定においても、差分法による地震動解析を有限要素法によるものに置き換え、幅広い問題設定に対する検討をより精度よく行えるようになることが期待されている。
内閣府政策統括官（防災担当）では、上記のように大規模有限要素法を用いた地震動計算の開発が先
行プロジェクトで進んでいたことに加えて、中央防災会議「東北地方太平洋沖地震を教訓とした地震・津波対策に関する専門調査会報告（H23.9）」において、揺れによる被害を軽減するための対策として、
長周期地震動対策等の必要性が指摘されたことや、「南海トラフ沿いの巨大地震による長周期地震動に関する報告（H27.12）」において、震源断層の極近傍に影響地域がある相模トラフは、長周期地震動の高層建築物等の影響評価に至らなかったことなどを政策的な背景として、「相模トラフ沿いの巨大地震に伴う長周期地震動による被害予測の高度化」の課題を進めてきた。その最終段階として、先行プロジェクトの大規模有限要素プログラムをベースとして本プロジェクトで整備した E-wave FEM を用いた長周期地震動の計算を「富岳」の政策枠で申請し、本年度計算を実施した。これは、国の被害想定のためのハザード評価の過程に先行プロジェクトからの開発の成果が社会実装されるという点で、本プロジェクトにとって重要な取り組みである。なお、内閣府防災は本プロジェクトの連携機関ともなっている。本プロ
ジェクトでは、内閣府防災・業務を担当する防災コンサルタント企業に対し、● プログラムの改良や使用方法の助言● E-wave FEM の計算精度等を示す計算結果の提供● 必要に応じた「富岳」の計算資源やデータ保存領域の拠出等を行い、政策枠を用いた計算をサポートした。②で述べる内容の一部は、この中の 2 点目に該当する。
その他、内閣府防災で実施された長周期地震動計算の精度を検証するため、本計算に対して 2 倍の分解能（データ量は約 8 倍）を設定したモデルでの計算を進めた。これにより、今後のプロジェクト遂行に必要となる、より高分解能な計算に必要な計算コストの情報も収集した。これは事業計画に記載した「弾性波動場の参照解」を作成したことに相当する。3 点目については、年度の終盤にはいって、計算結果の書き出しにおける不具合が政策枠での計算時に発生したことにより、本プロジェクトの計算資源やデー
タ保存領域を緊急避難的に提供することで、計算実施を支援することができた。
②E-wave FEM の精度検証と機能追加
本課題においては、①での大規模有限要素法での計算を進めるために必要な精度検証を行い、その結果を①での検討のために提供した。精度検証は、解析解との比較と、国の被害想定で用いられた実績のある既往手法との比較の双方により実施した。なお、本項目に関する数値計算には、富岳もしくは海洋研究開発機構の地球シミュレータを用いた。
②-(a) E-wave FEM の概要
E-wave FEM で用いられる計算手法について簡潔に述べる。E-wave FEM では弾性波動の支配方程式を有限要素定式化により離散化したシステムを扱う。すなわち
𝐊𝐮 + 𝐂𝐯 + 𝐌𝐚 = 𝐟-3-である。ここで𝐊、𝐂、𝐌、𝐮、𝐯、𝐚、𝐟はそれぞれ、剛性マトリクス、ダンピングマトリクス、質量マト
リクス、変位ベクトル、速度ベクトル、加速度ベクトル、外力ベクトルである。ダンピングマトリクス
は Rayleigh 減衰モデル[5]により設定される。時間ステップ幅に対する無条件安定性を重視し、陰的な時
間積分スキームである Newmark-𝛽法(𝛽=1/4, 𝛿=1/2)を適用することで、上記の方程式を𝐀𝛿𝐮𝑛 = 𝐛𝑛 ,
ただし、
𝐀=

4
2
𝐌 + 𝐂 + 𝐊,
2Δ𝑡Δ𝑡𝐛𝑛 = 𝐟𝑛 − 𝐪𝑛−1 + 𝐂𝐯𝑛−1 + 𝐌 (𝒂𝑛−1 +4𝐯 ),Δ𝑡 𝑛−1という形で数値的に時間積分する。ここで、δ𝐮、𝐪、𝐟はそれぞれ変位増分、内力、外力を表すベクトルである。Δ𝑡は時間増分であり、𝑛は時間ステップ数を示す。上記の連立一次方程式を解いた後、各ベクト
ル𝐪、𝐮、𝐯、𝐚はそれぞれ以下のように更新される。
𝐪𝑛 = 𝐪𝑛−1 + 𝐊𝛿𝐮𝑛 ,𝐮𝑛 = 𝐮𝑛−1 + 𝛿𝐮𝑛 ,𝐯𝑛 = 𝐯𝑛−1 +𝐚𝑛 = 𝐚𝑛−1 −

2
𝛿𝐮 ,
Δ𝑡 𝑛

4
4
𝐯
+
𝛿𝐮 .
Δ𝑡 𝑛−1 Δ𝑡 2 𝑛

毎時間ステップにおいて、巨大な疎行列𝐀（最も小規模なケースでも 107 オーダーの次元数）を含む線形方程式の求解が必要となる。共役勾配法・マルチグリッド法・精度混合演算を組み合わせた並列 solverを用い、線形方程式の求解を行う[3] [4]ことにより、「富岳」や京コンピュータなどのスーパーコンピュータを使用して多自由度（107-109 のオーダー、あるいはしばしばそれ以上）の問題を高速に解くことが可能となる。
精度検証の対象となる項目として、A.波動方程式の求解と B.震源モデリングに該当する精度検証を進めたので、その内容について述べる。表 1.2.1 に、精度検証・対応についての説明がある節を示す。
表 1.2.1 E-wave FEM の精度検証の概要
精度検証の対象とな導入されている手法

る事項
A.波動方程式の求解1.②において精度検証・対応についての説明がある節

有限要素法（非構造

速度波形：(d)現実的な地盤構造における長周四面体二次要素）・

期地震動計算の差分法との比較

Newmark-β 法

変位波形：(b)水平 2 層問題における変位計算の解析解との比較

B.震源断層のモデリ

点震源(c)面震源断層に対するスプリットノード法

ング

節点分割法[6]と多数点震源に基づく地震動計算結果の比較②-(b) 水平 2 層問題における変位計算の解析解との比較-4-

①の検討において、E-wave FEM で計算した変位波形を、地震計での観測データと比較する必要性が生じた。昨年度から進めてきている E-wave FEM の地震動計算の検証においては、速度波形のみを対象としてきた。そのため、本年度は水平二層問題における地震動計算の変位波形を E-wave FEM を用いて計算
し、波数積分法による解析解と比較することで、変位波形の計算精度を確認した。
堆積層が地表面から厚さ 2 ㎞、それ以深は母岩という設定を想定する。物性値、震源パラメータ、観測 点位置 を表 1.2.2、表 1.2.3、図 1.2.1 に 示す 。E-wave FEM で の計 算にお いては 、計算領 域を-40km<=x<=40km、-40km<=y<=40km、-54km<=z<=0km と設定する。2.5Hz 以下の周波数帯をターゲットとし、1 波長当たり最低 5 要素を確保する設定とした。地表面・層境界の水平方向メッシュサイズは 160mとした。解析解の計算には[9]の手法を用いた。このプログラムでは速度波形を出力するため、変位波形を数値積分により求めた。解析解・有限要素法の双方の計算結果を図 1.2.1 に示す観測点で比較した。こ
の設定下で、速度波形がよく一致することをあらかじめ確認した。
図 1.2.2 に計算結果を示す。E-wave FEM と解析解による変位波形は全体的によく合っているように見える。一方で、図中の拡大部分に示すように、地震波が通過した後の変位を見ると、両者に少しずれがあることがわかる。ここで、波形の一致度を簡便に評価するため、波形の相互相関係数をもとにしたス
コア𝐶 ∗を考える。すなわち、二つの時刻歴波形𝑎1 (𝑡)、𝑎2 (𝑡)を考えた時に𝐶 ∗ = 10max[𝐶(𝑎1(𝑡), 𝑎2 (𝑡)),0],
ここで𝐶(𝑎1 , 𝑎2 ) =∫ 𝑎1 (𝑡)𝑎2 (𝑡)𝑑𝑡2[∫ 𝑎1 (𝑡)𝑑𝑡]1/2 [∫ 𝑎22 (𝑡)𝑑𝑡]1/2である。𝐶 ∗ = 10であれば 2 つの波形は同一のものということになる。これは[10]の提案する波形のGOODNESS-OF-FIT を表す 10 のスコアのうち、
特に時刻歴波形の一致度に直接的にかかわるものとして、用いることとした。𝐶 ∗の計算時に 10 をかけているのは、GOODNESS-OF-FIT が波形の一致度に対して10 点満点の評価を与えることを目的として設定されていることに由来する（観測波形と計算波形の比較を想定した場合、4 未満は「poor」、4-6 は「fair」、6-8 は「good」、8 以上は「excellent」となるが、ここでは計算波形同士なので 10 に近い値での比較となる）。この指標をもとにすると、これらのずれは表 1.2.4の「No filter」の欄のように定量化される。このようなずれは長波長成分に出ていると考えられるため、カットオフ周波数 0.025Hz でハイパスフィルタをかけたところ、両者の変位波形のずれはかなり目立たなくなった(図 1.2.3) 。𝐶 ∗の指標によりずれの度合いを定量化すると、特に震源から遠方の点の z 成分など
で一致度の改善が顕著であることがわかる（表 1.2.4「Filtered」）。図 1.2.4「small model」に、本設定の計算結果について、最後の 5 秒間の変位の平均値を永久変位とみなして、各観測点での解析解との絶対差と相対差を比較したものを示す。ここで、絶対差:|E-wave FEM-解析解|, 相対差: |(E-wave FEM-解析解)/解
析解|と定義する。モデル端部に近づくにつれ、絶対差・相対差ともに大きくなることがわかった。そこ
で、FEM の水平計算領域を-80km<=x<=80km、-80km<=y<=80km と設定し、永久変位を同様に比較したところ、図 1.2.4 に示されるように、元の計算領域における計算結果に比べて両者の差が小さくなった。
これらのことから、E-wave FEM における計算領域端部の境界条件処理が、長波長での変位波形の計算
結果に影響していると推測される。計算領域を十分広くとることによりその影響を緩和することができるものの、速度波形がよく一致することを確認しただけの計算モデルでは、永久変位の計算精度が悪い場合があることが分かった。したがって、①で対象としている長周期地震動の計算で変位波形の結果を-5-用いる際には、計算領域を十分に広く確保するか、対象周波数の外側をフィルタで取り除いた波形のみを用いることが必要である。将来的に、変位波形の計算結果を、津波計算の初期値のような永久変位が必要な場合に活用するためには、こうした影響が出ないような境界条件の与え方を検討する必要があり、今後の課題である。
表 1.2.2 水平二層モデルの物性値
Thickness

Vp (m/s)Vs (m/s)

𝜌(km/m3)(km)
Sedimental

Quality factor

Quality factor

for P wave

for S wave

2

4,000

2,000

2,600

80

80

∞

6,000

3,464

2,700

140

140

layer
Bed rock

表 1.2.3 震源パラメータ
Hypocenter (km)Sourcetime

Rise time (s)Mo (N m)Strike, Dip, Rake2

1x10180o, 90o, 0ofunction shape(-3, 0, -6)

Triangular

図 1.2.1 観測点配置（黒点）と震央位置（赤バツ印）の関係-6-

図 1.2.2 フィルタなし変位波形の比較
左からそれぞれ、x, y, z 方向の変位波形を示す。凡例の「FEM」は有限要素解、「Hisada」は解析解、「disp=0」は比較用に記載した変位 0 に相当する線を表す。

表 1.2.4 抜粋した観測点における、変位波形各成分の C*スコアのフィルタ有り無しによる比較
Epicentral distance

6.40E+00

1.28E+01

1.92E+01

2.56E+01

3.20E+01x

1.00E+01

9.99E+00

9.94E+00

9.81E+00

9.59E+00y

1.00E+01

1.00E+01

9.97E+00

9.86E+00

9.62E+00z

1.00E+01

1.00E+01

9.78E+00

9.48E+00

8.90E+00x

1.00E+01

9.99E+00

9.96E+00

9.85E+00

9.66E+00y

1.00E+01

1.00E+01

9.98E+00

9.94E+00

9.83E+00z

1.00E+01

9.99E+00

9.96E+00

9.96E+00

9.94E+00(km)

No filter

Filtered-7-図 1.2.3 カットオフ周波数 0.025Hz のハイパスフィルタをかけた変位波形の比較
左からそれぞれ、x, y, z 方向の変位波形を示す。凡例の「FEM」は有限要素解、「Hisada」は解析解を表す。
「disp=0」は比較用に記載した変位 0 に相当する線を表す。

-8-

図 1.2.4 フィルタをかけない変位波形について、最後の 5 秒間の変位の平均値を永久変位とみなして、各観測点での両者の絶対差（|E-wave FEM-解析解|）と相対差（|(E-wave FEM-解析解)/解析解|）を比較したものを示す。Small model は図 1.2.1Large model は水平計算領域を広げた結果に基づくもの。ｚ成分の相対差については、そもそもの応答レベルが他成分と比べて小さいために相対差が比較的大きいこと、また応答の符号が途中で入れ替わることで局所的に相対差が大きくなる部分あり、それが縦軸の最大値を超えていることに注意する必要がある。

②-(c) 面震源断層に対するスプリットノード法と多数点震源に基づく地震動計算結果の比較①では、面的な広がりをもった震源モデルに対する応答に対する計算が必要となる。面的な広がりをもった震源モデルを差分法、有限要素法などの数値計算に導入する方法としては、(i) 一定の間隔で分布する点震源の集合体を一つの震源モデルとみなす方法（多重点震源）(ii) 節点分割法[6]の二つが考えられる。節点分割法は、断層面上に位置する節点の両側での仮想的な食い違い変位を考え、それと等価な体積力に対する応答を計算する方法である。有限要素法において、非構造要素により断層
面位置節点を配置して節点分割法を導入することで、点震源の集合を考える場合と比べて、断層面上の連続的なすべりをより精度よく考慮することができるものと期待されている。そのため①においても、
節点分割法に基づいた震源モデルの導入を前提とした検討が進められている。E-wave FEM では、(i)(ii)の双方を扱うことが可能となっており、両者を比較することで、震源断層のモデリングに関する精度検証を行う。昨年度も類似の比較を行い、本課題の令和二年度成果報告書 2-2(1)①-c「面的な広がりを持った震源モデルに対する地震動計算の精度検証」にその結果について記載した。ここでは、断層の両側で物性値の変わらない平面断層での面震源に対して(i)と(ii)の比較を行い、両者の計算波形に明確な差がみ-9-られた。しかし本年度の調査で、記載した結果を計算した際の節点分割法の計算設定ファイルに誤りがあることが分かった。この誤りにより、面震源上の破壊伝播方向の x・y 成分が想定したものと逆になった状態で計算した結果を記載していた。そこで本年度は、まず設定ファイルの誤りを正し、面震源上の破壊伝播方向が正しく設定されるようにした。その上で、①で必要とされる、より一般的な設定下での検討として、「面震源が曲面に入力される場合」と「面震源が物性（速度）境界に入力される場合」の二つの設定で検証を行った。
「面震源が曲面に入力される場合」について、半無限弾性体中に設置した曲面に面震源を入力する設定を考える。物性値、震源パラメータ、観測点・断層面位置を表 1.2.5、表 1.2.6、図 1.2.5 に示す。曲面
は、xy 平面上で正方形、xz 平面上でサインカーブ型（深さ 6km 程度）をしている。すべり量は x・z 成
分のみを持ち、面に沿う方向に下盤側から見て一様に大きさ 1ｍとする。図中星印に示される破壊開始点
から 3km/s で同心円状に破壊伝播する。
E-wave FEM での計算においては、計算領域を-22.5km<=x<=22.5km、-22.5km<=y<=22.5km、-17km<=z<=0km と設定する。2.5Hz 以下の周波数をターゲットとし、1 波長当たり
最低 5 要素を確保する設定とした。地表面・層境界の水平方向メッシュサイズは 125m とした。多重点震源計算については、曲面に沿って等間隔になるように 2,548 個の点震源を 125m 間隔で配置した。双方の計算結果を図 1.2.5 に示す観測点で比較した。図 1.2.6 に速度波形の比較を示す。節点分割法と多重点震源による計算結果は全体的によく合っているように見える。一方で、図中の中央拡大部分に示すように、波形を拡大すると、振幅にわずかなずれがある部分があることが分かった。そこで、地表面・層境界の水平方向メッシュサイズを 50m とし、多重点震源計算を 15,488 個おいて比較したのが図 1.2.7 である。
ここでは、メッシュサイズが大きい場合に比べ、振幅のずれが改善されている。ただし、メッシュサイ
ズの変更による波形一致度の改善を𝐶 ∗スコアにより定量化すると、x 成分のいくつかの観測点ではメッシ
ュを細かくしたことにより一致度が悪くなっている（表 1.2.7）。メッシュサイズ変更により曲面のモデル化がわずかに異なっていることによる影響を比較の中で切り分けきれていないことが原因である可能性がある。もっとも、メッシュサイズの変更前後双方で波形全体の一致は非常に良いといえ、今回の検証目的に対して十分な検討ができていると考える。
「面震源が物性（速度）境界に入力される場合」について、水平二層構造の速度境界を断層とみなした設定を考える。物性値を表 1.2.8 に示す。震源パラメータと観測点位置は表 1.2.5、図 1.2.5 上と同様とした。断層は速度境界に位置し、xy 平面上で正方形（図 1.2.5 上と同様）の平面断層（深さ 6km）である。
すべり量は x 成分のみを持ち、面に沿う方向に下盤側から見て一様に大きさ 1ｍとする。図中星印に示される破壊開始点から 3km/s で同心円状に破壊伝播する。E-wave FEM での計算においては、計算領域を-24km<=x<=24km、-24km<=y<=24km、-27km<=z<=0km と設定する。2.5Hz 以下の周波数をターゲットとし、1 波長当たり最低 5 要素を確保する設定とした。地表面・層境界の水平方向メッシュサイズは 75mとした。多重点震源計算については、曲面上で等間隔になるように 6,561 個の点震源を 75m 間隔で配置した。面震源におけるすべり量を各点震源の地震モーメントに変換するためには剛性率が必要であるが、面震源が物性境界に位置する場合においては境界の両側のいずれの物性値も直接用いることができない。
本検討では、[11]により示される「有効剛性率」、すなわち両側の剛性率の調和平均を用いることとした。
図 1.2.8 に速度波形の比較を示す。節点分割法と多重点震源による計算結果は全体的によく合っているよ
うに見える。一方で本設定においても、図中の中央拡大部分に示されるように、波形を拡大すると、振幅にわずかなずれがある部分があることが分かった。そこで、地表面・層境界の水平方向メッシュサイ- 10 -

ズを 25m とし、多重点震源計算を 58,081 個おいて比較したのが図 1.2.9 である。ここでは、波形を拡大しても振幅がよく一致している。メッシュサイズの変更による波形一致度の改善は、𝐶 ∗ スコアにより定量化した（表 1.2.9）。
「面震源が曲面に入力される場合」と異なり、すべての成分・観測点で一致度が改善している。本設定においては断層が平面であり、メッシュサイズによらず断層のジオメトリを完全にモデル化できていることによる違いである可能性がある。ただここでも、メッシュサイズの変更前後双方で波形全体の一致は非常に良いといえ、いずれの場合にしても今回の検証目的に対して十分な検討ができていると考える。
これより、「面震源が曲面に入力される場合」と「面震源が物性（速度）境界に入力される場合」の双方において、適切な設定下では節点分割法と多重点震源による計算結果はよく一致することが分かった。
一方で、点震源の配置の密度が十分でない場合には両者に振幅のわずかなずれが発生し得ることも確認された。従来広く使われる多重点震源の代わりに節点分割法を用いることの利点をより詳しく整理する
ことは今後の課題である。
表 1.2.5 半無限弾性体モデルの物性値
Vp (m/s)

6,000

𝜌(km/m3)

Vs (m/s)

3,464

2,700

Quality factor

Quality factor

for P wave

for S wave

9,999

9,999

表 1.2.6 震源パラメータ
Source time function shape

Rise time (s)Triangular2- 11 -

図 1.2.5 半無限弾性体における観測点配置（黒点）、断層面（赤四角・線）、破壊開始点（青バツ印）の関係- 12 -

図 1.2.6 半無限弾性体での曲面断層に対する速度波形の比較（水平方向メッシュサイズ 125m）
左からそれぞれ、x, y, z 方向の速度波形を示す。

図 1.2.7 半無限弾性体での曲面断層に対する速度波形の比較（水平方向メッシュサイズ 50m）
左からそれぞれ、x, y, z 方向の速度波形を示す。

- 13 -

表 1.2.7 半無限弾性体での曲面断層に対する、抜粋した観測点における速度波形各成分の C*スコアの、
地表面・層境界の水平方向メッシュサイズ変更の前後の比較
Epicentral distance (km)2.9699E+00

5.9397E+00

8.9096E+00

1.1879E+01

1.4849E+01

x9.9955E+00

9.9974E+00

9.9990E+00

9.9984E+00

9.9970E+00

y9.9928E+00

9.9990E+00

9.9994E+00

9.9990E+00

9.9969E+00

z9.9987E+00

9.9988E+00

9.9981E+00

9.9987E+00

9.9990E+00

x9.9991E+00

9.9953E+00

9.9949E+00

9.9975E+00

9.9986E+00

y9.9980E+00

9.9992E+00

9.9998E+00

9.9993E+00

9.9996E+00

z9.9989E+00

9.9996E+00

9.9997E+00

9.9998E+00

9.9998E+00

Quality factor

Quality factor

for P wave

for S wave

125m

50m

表 1.2.8 水平二層モデルの物性値
Thickness

Vp (m/s)Vs (m/s)

𝜌(km/m3)(km)
Sedimental

6

4,000

2,000

2,600

80

80

∞

6,000

3,464

2,700

140

140

layer
Bed rock

図 1.2.8 速度境界面震源に対する速度波形の比較（水平方向メッシュサイズ 75m）
左からそれぞれ、x, y, z 方向の速度波形を示す。

- 14 -

図 1.2.9 速度境界面震源に対する速度波形の比較（水平方向メッシュサイズ 25m）
左からそれぞれ、x, y, z 方向の速度波形を示す。

表 1.2.9 速度境界面震源に対する、抜粋した観測点における速度波形各成分の C*スコアの、地表面・層境界の水平方向メッシュサイズ変更の前後の比較
Epicentral distance

2.9699E+00

5.9397E+00

8.9096E+00

1.1879E+01

1.4849E+01

X

9.9988E+00

9.9955E+00

9.9920E+00

9.9945E+00

9.9975E+00

Y

9.9980E+00

9.9983E+00

9.9882E+00

9.9937E+00

9.9978E+00

Z

9.9986E+00

9.9971E+00

9.9968E+00

9.9991E+00

9.9977E+00

X

9.9999E+00

9.9996E+00

9.9994E+00

9.9996E+00

9.9998E+00

Y

9.9998E+00

9.9999E+00

9.9990E+00

9.9995E+00

9.9998E+00

Z

9.9998E+00

9.9997E+00

9.9997E+00

9.9999E+00

9.9998E+00

(km)

75m

25m

②-(d) 現実的な地盤構造における長周期地震動計算の差分法との比較①の最初に述べたように、国の被害想定のための長周期地震動計算には、差分法による地震動計算手法が用いられてきた。これを E-wave FEM による計算に置き換えることを検討するにあたり、現実的な地盤構造に基づく問題設定の下で両者を比較することは重要である。国の被害想定のための地震動計算においては、差分法での計算において多く適用された実績のある、地下の現実な地盤構造は考慮するが地形は考慮せず平坦とした設定が採用されてきた。このような設定であれば、「有限要素法と差分法の結果が一致すべき」である。このような設定下で両者の数値解を比較し、差異が合理的な範囲にとどまって- 15 -

いるかを確認する。
まず強震動評価に用いる差分法による三次元地震波動伝播解析コードを用いて、現実的な地盤構造を対象とした地震動計算を行った。地震波動伝播解析は、運動方程式と等方弾性体を仮定した構成方程式（フックの法則）を差分方程式に変換し、弾性波動方程式の速度と応力の時間発展を数値的に解くことにより行った（例えば[12]）。空間 4 次、時間 2 次の差分近似精度とした食い違い格子を用いている。非弾性減衰の効果は、[13]の方法に基づいて考慮した。振幅の減衰性を考慮する参照周波数を f0、時間ステップを Δt、クオリティ・ファクターを Q(x,y,z)とすると、各時間ステップの速度と応力を計算する際に、次式で表される係数 A(x,y,z)を 1 つ前の時間ステップの値に乗じることで簡易に非弾性減衰の効果が評価される。

  f 0 t Ax, y, z   exp
 Q  x, y , z  （1.2.1）

この方法では P 波の減衰 Qp と S 波の減衰 Qs を独立に与えることはできない。また、Q 値の周波数依存性として Q  Q0  f が仮定される。計算領域の側面と下面では、有限の領域を打ち切るために、[14]のone-way 波動方程式に基づく無反射境界条件と、[8]の吸収境界条件を併用した境界条件を設けた。[8]による吸収境界は、吸収率 W を W  exp   2  J 0  j 2 で表し、境界に近づくに従い指数的に W を増加させる。α は吸収率の大きさをコントロールする係数、J0 は吸収領域幅を表す。J0 を大きくするほど効果は大きく、対象とする計算問題に応じて J0 の幅を決める。計算領域の上面では、[15]による Improved Vacuum
Condition を採用し、任意形状の自由地表面を考慮できるようにしている。[16]の方法に基づき、応力の各成分に対して、対応するモーメントテンソルを入力することにより、食い違い震源断層を考慮する。
計算効率の向上のために、領域分割に基づいた並列計算を行った。3 次元の計算領域を NS、EW、UDの 3 方向に分割し、分割した小領域の演算を複数のノードに割り当てて計算負荷を分散させた。計算ス
テップ毎に MPI を用いた非同期通信を行い、各小領域の袖領域に設けたオーバーラップグリッドの速度と応力の値を、隣接するプロセス間で送受信して共有した（例えば、[17][18]）。ノード内の並列計算もプ
ロセス並列とした flat-MPI により行う。
計算モデルは、[19]に準拠して設定した。[19]では、関東平野について、2009 年伊豆半島沖の地震（M5.1）を対象とした地震動シミュレーションを行っている。地盤構造の不整形性に加えて、地形、さらに、海水を考慮することで、計算される地震動特性に違いが生じる結果を示している。ここでは、地形と海水を考慮せず、地盤構造の不整形性のみを考慮した地震動計算を行った。計算に用いた地震本部[20]の地盤モデルについて、物性値を表 1.2.10 に、地震基盤（S 波速度 Vs＝3.2 km/s 層）の上面深度分布を図 1.2.10に示す。図中の●印の位置で、2009 年伊豆半島沖の地震（M5.1）の観測記録が得られている。水平方向の計算領域は、NS（X）方向 190 km、EW（Y）方向 140 km の範囲（図 1.2.10 の図示範囲）とした。深さ（Z）方向では、深さ約 60 km までを計算領域とした。格子間隔は、強震動評価や地震動シミュレーシ
ョンに際して、一般的に用いる計算条件として、水平 0.1 km、上下 0.05 km（ただし、深さ 10 km 以深は0.2 km）とした。計算領域の東西南北の端部および下端に、240 grid 幅の吸収領域を設けた。格子点数は（X, Y, Z）=（2380, 1880, 736）grid である（計約 33 億）。計算は、地球シミュレータ ES4 を用いて、X、Y、Z の 3 方向をそれぞれ 10、8、8 分割した 640 並列により行った。格子間隔が計算される地震動に与える影響について検討するために、格子間隔を水平 0.05 km、上下 0.02 km（ただし、深さ 10km 以深は- 16 -

0.05km）とした地震動計算を別途行った。計算領域の東西南北の端部および下端に、480grid 幅の吸収領域を設けた。格子点数は（X, Y, Z）=（4760, 3760, 2184）grid である（計約 391 億）。計算は、地球シミュ
レータ ES4 を用いて、X、Y、Z の 3 方向をそれぞれ 35、16、8 分割した 4,480 並列により行った。時間
刻みは共通して 0.001 秒、時間長は 200 秒とした。震源として、図 1.2.10 の☆印の位置に、図 1.2.11 に示
す F-net のメカニズム情報に基づく点震源を与えた。震源時間関数は、図 1.2.12 に示すライズタイムを 0.9秒とした二等辺三角形型で与えた。観測記録が得られている地点（図 1.2.10●印）に隣接する格子点で計算結果を出力した。また、海底で計算される地震動特性を調べるために、図 1.2.10▲印の位置に出力点を設けた。計算結果の上下成分は、UP を+、Down を－となるように符号を変えた。また、計算結果に対して、図 1.2.13 に示す周期特性のバンドパスフィルター処理を行った。
格子間隔を水平 0.1km、上下 0.05km とした計算結果として、震央と東京湾岸をつなぐ図 1.2.14 の測線について、最大振幅で基準化した NS 成分、EW 成分、UD 成分の速度ペーストアップを図 1.2.15 に示す。
×印と●印は、JMA2001 速度構造[21]に基づく P 波と S 波の到達時刻を表す。格子間隔を水平 0.05 km、上下 0.02 km に狭めた計算結果として、図 1.2.15 と同様の速度ペーストアップを図 1.2.16～18 に示す。dx、dy、dz は、NS 方向、EW 方向、上下方向の格子間隔を表す。格子間隔を狭めたことにより、計算結果に違いが生じる。
次に、E-wave FEM を用いて現実的な地盤構造における地震動計算を行い、差分法（FDM）を用いた地震動計算結果との比較を行った。モデルは FDM を用いた計算と同様に[19]に準拠して設定した。検討対
象領域を図 1.2.10 に示した東西 190km×南北 140km の領域に東西南北それぞれに 24km 幅の吸収帯を加えた東西 238km×南北 188 km の領域とし、深さ方向は、48km 幅の吸収帯を含み 108km とした。観測点数は226 点である。地盤モデルの設定パラメータは表 1.2.11 のものを用いた。メッシュについては、最小メッシュサイズ DS とメッシュ数 NX、NY を指定する形とした。FMAX は最小 Vs の 350 m/s および DS から
0.7 Hz とした。物性値のうち Q 値については Qs の値を用いることとした。図 1.2.19 に可視化した第 14、18、21 層上端の形状を示す。総要素数は約 17.7 億、総節点数は約 23.8 億となった。
計算時間は 0.01 秒刻みで 200 秒間とした。
総ステップ数は 20,000 である。
震源は FDM と同様に図 1.2.10に☆印で示す位置の深さ 5km 地点に、F-net を参照して、M0＝2.81×1016Nm、（strike、dip、rake）＝（168°、
90°、12°）の点震源を置いた。震源時間関数も同様にライズタイム 0.9 秒の 2 等辺三角形とした。波形には図 1.2.13 に示す周期 2~10s のバンドバンドパスフィルターを施した。
計算はプリプロセスに海洋研究開発機構の Jasper、本計算およびポストプロセスに「富岳」を用いた。本計算は 385 ノード 1,540MPI プロ
セス数で行った。プリプロセス、本計算、ポストプロセスの各ステップの実行にかかった計算時間を表
1.2.12 に示す。E-wave FEM と FDM の比較を行った観測点は図 1.2.15 に示されたものである。FDM の計
算は格子間隔の異なる 2 通りの計算を実施しているため、それぞれについて FEM との比較を行った。
一般的な格子間隔の場合（水平 100m 上下 50m）の FDM との、計算速度波形の比較を図 1.2.20 に、フーリエスペクトルの比較を図 1.2.21 に示す。また、格子間隔を細かくした場合（水平 50m 上下 20m）のFDM との、計算速度波形の比較を図 1.2.22 に、フーリエスペクトルの比較を図 1.2.23 に示す。バンドパ
スフィルターをかける際は、震源時（図 1.2.22 の 0 秒）から 200 秒間を切り出し、0.3 秒の cos テーパーを最初と最後に施した。速度フーリエスペクトルの計算も同様に 200 秒間の波形を用いた。まず、一般的な格子間隔の場合の FDM と FEM とを比較する。速度波形（図 1.2.20）を見ると、実体波と見られる- 17 -

部分は震央距離 50km 程度まででは 2 手法で類似した結果が得られているものの、それより遠い距離では位相は似ているものの振幅に差が見られる。また、表面波とみられる部分については、最も近い観測点
AAA001 においても分散性や振幅に違いが見られ、すべての観測点で実体波と思われる部分に比べて差が大きい。フーリエスペクトル（図 1.2.21）を見ると、各観測点で類似した形状となっているが、スペクトルのピーク周波数に違いが見られる他、振幅についても特に震源から遠い観測点で差が比較的大きい。
また、FEM のほうが FDM に比べスペクトル振幅が大きい傾向にある。
格子間隔を細かくした FDM と FEM を比較すると、速度波形（図 1.2.22）を見ると、一般的な格子間隔の場合（図 1.2.22）に比べ、表面波と見られる部分の位相のずれがかなり改善されている。フーリエス
ペクトル（図 1.2.23）においては、FEM のほうが FDM に比べスペクトル振幅が大きい傾向にあるのは一般的な格子間隔の場合と同様だが、スペクトルの山谷の周波数が両手法で一致している。本検討に用いた地盤モデルは最小 Vs が 350m/s であり、対象周波数が 0.5Hz までである。そのため、一般的に用いられる 1 波長 5 グリッド／メッシュの条件を考えると、必要なグリッド／メッシュサイズは 140m となり、本検討の FEM と 2 つの FDM 計算のメッシュ／グリッドサイズは十分小さい。しかし、2 つの FDM の結果を、FEM を通して比較すると、表面波と思われる部分で違いがあり、フーリエスペクトルのピーク周波数にも差が見られた。おそらく、メッシュサイズの違いによる主に表層部分のモデル化の違いの影響と考えられる。FEM とメッシュが細かい場合の FDM においてもまだ差異はあるが、その要因としては、
層境界のモデル化の違いや減衰の設定の違いが考えられる。FEM と FDM の結果をある程度一致させる
ためには、FDM のメッシュ間隔を一般的な間隔より細かくする必要があると考えられるが、その程度については許容する差異によるため、定量的な検討が今後必要であると考えられる。
ここでは、E-wave FEM および差分法を用いて、関東平野の不整形な地下構造と伊豆半島沖の点震源に対する、震央から東京湾岸をつなぐ測線の地表面 15 地点の地震動を計算し、速度波形、フーリエスペク
トルについて両手法の結果を比較した。一般的な格子間隔の場合と比べて、格子間隔を細かくした場合の差分法の結果は、表面波と見られる部分の E-wave FEM との位相のずれがかなり改善され、フーリエス
ペクトルの山谷の周波数も E-wave FEM と一致した。これより、従来使われる地震動数値計算の手法と、E-wave FEM の計算結果は、一致すべき問題設定下において適切な計算条件を用いれば、よく一致することが示されたといえる。2 つの差分法の結果の違いは、格子間隔の違いによる主に表層部分のモデル化の違いの影響と考えられる。E-wave FEM と格子間隔が細かい場合の差分法においてもまだ差異はあるが、その要因としては、層境界のモデル化の違いや減衰の設定の違いが考えられる。E-wave FEM と差分法の結果をある程度一致させるためには、差分法の格子間隔を一般的な間隔より細かくする必要があると考えられるが、その程度については許容する差異によるため、定量的な検討が今後必要であると考えられる。それに関連して、格子間隔をさらに細かくした差分法の検討や、プレポストの時間を含めた E-wave
FEM と差分法の性能の比較も必要であると考えられる。加えて、層境界のモデル化の違いを可視化して把握することも重要だと考えられる。また、本検討では①に合わせて地下構造の不整形性のみを考慮したが、それに加えて地形を考慮した場合や地形と海水を考慮した場合の検討、さらに震源についても別の点震源や面震源に対する検討が、現実的な地盤構造への長周期地震動計算手法の適用性を確認する上で重要と考えられる。

- 18 -

表 1.2.10 地盤モデルの物性値- 19 -

●印：K-NET、 KiK-net、 MeSO-net 観測点▲印：海底の計算結果確認のために設定した出力点☆印：震央
図 1.2.10 S 波速度 3.2 km/s 層の上面深度分布- 20 -

図 1.2.11 F-net によるメカニズム情報

図 1.2.12 震源時間関数（振幅基準化）

図 1.2.13 バンドパスフィルターの周期特性- 21 -

図 1.2.14 波形ペーストアップに用いた観測点（○印）
赤色のコンターは、S 波速度 3.2 km/s 層の上面深度を表す- 22 -

図 1.2.15 周期 2～10 秒の速度ペーストアップ（振幅基準化）
×印と●印は、JMA2001 速度構造に基づく P 波と S 波の到達時刻を表す- 23 -

格子サイズ：水平方向 100 m、 鉛直方向 50 m

格子サイズ：水平方向 50 m、 鉛直方向 20m

図 1.2.16 NS 成分の速度ペーストアップの比較（周期 2～10 秒、振幅基準化）
×印と●印は、JMA2001 速度構造に基づく P 波と S 波の到達時刻を表す- 24 -

格子サイズ：水平方向 100 m、 鉛直方向 50 m

格子サイズ：水平方向 50 m、 鉛直方向 20m

図 1.2.17 EW 成分の速度ペーストアップの比較（周期 2～10 秒、振幅基準化）
×印と●印は、JMA2001 速度構造に基づく P 波と S 波の到達時刻を表す- 25 -

格子サイズ：水平方向 100 m、 鉛直方向 50 m

格子サイズ：水平方向 50 m、 鉛直方向 20m

図 1.2.18 UD 成分の速度ペーストアップの比較（周期 2～10 秒、振幅基準化）
×印と●印は、JMA2001 速度構造に基づく P 波と S 波の到達時刻を表す- 26 -

表 1.2.11 地盤モデルの設定パラメータ
項目

説明

値

NS

地盤層の数

23

NK

各地盤層における有限要素分割の階層の深さ

3SF

スクオッシュファクター

2

DS

メッシュの水平方向最小サイズ

100 m

NX

x 方向のメッシュ数

1880

NY

y 方向のメッシュ数

2380

NOW

DS で指定されたメッシュサイズに含まれる波数

5

FMAX

離散化及び Rayleigh 減衰の対象とする最大周波数

0.7 Hz

図 1.2.19 可視化した第 14、 18、 21 層の上面形状
表 1.2.12 各ステップの計算時間
Step

総計算時間 (h)プリプロセス（メッシュ生成・領域分割）

33.2

地震動計算

7.6

ポストプロセス（可視化用データ生成）

2.8- 27 -

図 1.2.20 FEM と一般的な格子間隔の FDM の計算速度波形の比較（各波形ペアで振幅を基準化）- 28 -

NS 成分、EW 成分、UD 成分

図 1.2.21 FEM と一般的な格子間隔の FDM のフーリエスペクトルの比較- 29 -

図 1.2.22 FEM と格子間隔を細かくした FDM の計算速度波形の比較（各波形ペアで振幅を基準化）- 30 -

NS 成分、EW 成分、UD 成分

図 1.2.23 FEM と格子間隔を細かくした FDM のフーリエスペクトルの比較- 31 -

③E-cycle FEM を用いて計算した弾性応答グリーン関数の数値誤差定量化
地殻変動計算プログラム E-cycle FEM を用いて計算できるプレート境界における断層すべりに対する
静弾性グリーン関数は、地震時の断層すべり分布推定や地震間の断層固着度推定に加え、津波シナリオ
計算の際の海底地殻変動推定など、ハザード評価に資する検討に対して様々の用途がある。そのため、沈み込み帯地震に関する地震学・地震工学研究における重要なコンポーネントと位置付けられる。線形性により、一度計算しておけば重ね合わせによりいろいろな現象の解析に適用できる。半無限弾性体のような単純な媒質を仮定する場合においては[22][23]に代表される計算プログラムの公開により、弾性グ
リーン関数計算が広く行われ、測地観測と組み合わせた断層すべり逆解析などの分野で成果を挙げてきた。一方、プレート境界周辺の媒質に関しては、地表面形状はもとより、地下の弾性構造についても、
地震探査やトモグラフィにより情報が蓄積されるようになってきた。それらの知見に基づき、地下の地
震波速度構造のデータベースが整備されてきている（例えば[24][25][26][27]）。また地上の GNSS 地殻変動観測が常時行われている(例えば the GNSS Earth observation network system (GEONET)[28])ことに加え、海底地殻変動観測（GNSS-acoustic (GNSS-A)観測[29][30][31][32][33]・海底水圧計[34][35][36][37])や干渉
合成開口レーダー (InSAR)[38][39][40]による地殻変動解析がより多くの地域で可能となってきた。これらの観測から得られるデータがカバーするプレート境界周辺の広い範囲での地殻変動データを整合的に扱おうとする場合、地形、弾性パラメータ（例えば[41][42][43][44][45][46][47][48]）、地球の楕円性[49]を考慮したグリーン関数の計算、および断層すべり逆解析、津波計算などが望ましいと考えられるようになってきた。計算結果の信頼性の観点から、公表されている地下の弾性構造データベースや、標準として用いられている測地系に対して可能な限り忠実（high ﬁdelity）であることも重要である[50]。このような計算には、任意の形状を取り入れやすく、地表面での応力フリーの境界条件を自動的に満たす有限要素法が向いているため、有限要素法を用いた断層すべりに対する静弾性グリーン関数の計算例がここ 10
年くらいで増えてきている。しかし、有限要素法による計算は計算モデル生成・弾性応答計算ともに、解析的な表現が可能な半無限弾性体における応答計算と比べ、はるかに計算コストが大きい。オペレー
ションについても、半無限弾性体での計算に比較すると、ある程度の訓練が必要となる。また、グリーン関数には、一度計算しておけば後はユーザが簡便な重ね合わせ計算によりいろいろな現象の解析に活用できるという特徴がある。そのため、一度計算されたグリーン関数は、ライブラリとして計算者とは異なるユーザに様々な用途に使われることが往々にしてあると考えられる。たとえば、[50]によりライブ
ラリとして発表されている、南海トラフ域と日本海溝域を対象とした、3 次元弾性不均質性を考慮した有限要素計算によるグリーン関数はその一例である。公開されるされないにかかわらず、今後も、国内外の主要な沈み込み帯・地震発生帯において同様のグリーン関数ライブラリが計算されていくものと思われる。
数値シミュレーションにおける計算パラメータの設定に任意性がある中で、グリーン関数ライブラリのような科学的・工学的な成果物における数値解の信頼性を担保することは重要である。数値解の信頼性を評価する枠組みとして Veriﬁcation & Validation という概念が提案されている。Veriﬁcation とは、ある物理現象に対して想定される数値モデル（偏微分方程式など）が、シミュレーションコードによって正しく解かれているかどうかを確認する作業である。Veriﬁcation は計算コード、計算アルゴリズム、計算モ
デル（計算設定）、誤差の定量化までの確認を含んでいる。Validation は、数値モデルが対象とする物理現象を合理的に説明できるかどうかを確認するプロセスである。このように分類することで、数値解の結- 32 -

果と実際の観測データとの差異を議論する際の混乱を回避することができる。工学分野では、その仕様標準までが整備されている[51]。一方で、グリーン関数ライブラリの主な適用先である地震学の分野での学術的研究においては、V＆V の概念は浸透していない。例外はハザード評価につながる Engineering
Seismology の分野における veriﬁcation[52]や、広く使われるシミュレーション手法に対して学術コミュニ
ティが中心となって進める veriﬁcation[53][54]などである。ただ、ここでの主な目的は計算コードやアルゴリズムの検証と手法間の比較であり、数値誤差の定量化まで至っている研究はほとんどないのが現状である。
一方で、グリーン関数ライブラリはユーザ各自による重ね合わせ計算により二次利用されることが期待されるので、ライブラリを作る段階でどのような用途に使われるかがわからない。また、グリーン関数を計算した人が自分でそれを解析に用いる場合も、同じグリーン関数をライブラリと同様の形で様々な現象に適用することは多いはずである。そのため、グリーン関数ライブラリに関する使用限界等について、使用する際にあらかじめ判断できる材料があることが望まれる。その指標とするための情報として、数値計算に含まれる誤差をあらかじめ定量化しておくことが望ましい。多くの場合、グリーン関数
ライブラリについて veriﬁcation を行う対象となるのは、メッシュサイズが十分小さいかということと、地球に相当する楕円体から計算領域を切り出すことによる、領域端部での境界条件の影響が挙げられる。
[50]はグリーン関数ライブラリの計算に際して、過去の類似の計算設定における検討例[55]を参考にしつつ、選択したいくつかのケースに対して、モデルパラメータの変化に対して計算結果が大きく変わらないことを数回の予備的解析によって確かめたが、その内容を論文の中で明示せず、定量的な議論も行っていない。このグリーン関数ライブラリの作成においては、50 億自由度程度の計算モデルを用いて、南海トラフ域と日本海溝域それぞれについて、1,000 回弱の静弾性計算が実施された。これだけでも、例えば前者については JCAHPC[1]の Oakforest-PACS[57]（1 ノード辺り理論性能 3TFLOPS 程度）を 512 ノードを用いて 16 時間余りの計算が必要であるように、計算コストはかなり大きいものである。これに比較し、メッシュをさらに細かくしたり、領域サイズ をさらに大きくしたりした計算を、1,000 回弱の全計算ケースに対して網羅的に計算して誤差評価をすることは、計算時点ではあまり現実的でないと考えていた。このように、ターゲットとする設定よりさらに大規模な計算モデルを参照した veriﬁcation が、計算コストの問題でしばしば難しくなることは、[58]でも指摘されている。しかし、「富岳」の登場によって、より細かいメッシュや大きな領域での計算結果の veriﬁcation や数値誤差の定量的評価を系統的に行うことが現実的になった。[50]のグリーン関数ライブラリだけでなく、数値計算ベースのグリーン関数が広く計算され、利用されることを考えると、このようなより大規模な計算環境の使用を前提としたグリーン関数の数値誤差評価手法を検討することは意義深いと考える。
これらの背景に基づき、本項目では、有限要素法等の数値計算手法に基づいて計算されるグリーン関
数に対し、「富岳」という大規模計算環境を活用し、よりメッシュを細かくした・領域を広くした計算モデルを用いて、数値誤差の定量化を[50]のグリーン関数ライブラリに対して行った。数値誤差の定量化の結果と、グリーン関数ライブラリが適用されることが想定される実現象の特性とを比較することで、計算精度の観点からのグリーン関数ライブラリの適用性の限界を議論し、利用者のアプリケーション選択に役立てる。このように、最先端の大規模計算機とアルゴリズムを活用して、より大規模な計算モデルの数値解を参照することで、高詳細モデルを用いて計算したグリーン関数ライブラリの検証と数値誤差の定量化を試みた。この「より大規模な計算モデルの数値解」の計算により、事業計画に記載した「弾- 33 -

性地殻変動応答の参照解」を作成したことになる。
「グリーン関数ライブラリ」において扱うグリーン関数とは、対象とするプレート境界を離散化する
小断層に対する弾性応答変位である。小断層に対する応答変位を計算する数値解析手法として、地殻構
造のような形状の複雑な連続体の解析に適した有限要素法を用いる。地殻変動を静弾性変形としてモデル化する。境界条件を考慮して整理すると、問題は以下の一次連立方程式に帰着する。
𝐊𝐮 = 𝐟
ここで K、u、f はそれぞれ全体剛性マトリクス、変位ベクトル、外力ベクトルである。小断層におけるすべり量を節点分割法[1]により外力ベクトル f に変換する。計算モデルの側面と底面にはディリクレ
条件が適用される。具体的には、表面の法線ベクトルの方向の変位を 0 とし、それ以外の方向には変位を与えない、という条件である。楕円体や半無限領域から対象領域を切り出し、その領域端部に変位固定の境界条件を課している。これらのモデル端部での条件が計算結果にどの程度影響を及ぼしているかを定量的に把握する必要がある。
本式は巨大な線形連立方程式であり、代表的な反復求解アルゴリズムの一つである共役勾配法により求解する。求解の収束性を、四面体メッシュに基づく一次要素を用いた幾何的マルチグリッド法により、スケーラビリティを保ちつつ大幅に改善している[44]。ここに代数的マルチグリッド法も導入することで収束性がさらに改善され、1012 オーダーの自由度の計算モデルに対しても、現実的な時間での求解が可能となる[60]。本手法において可能な計算の詳細さは、我々の知る限り、固体地球分野において他に類を見ない。
有限要素モデルは、[44][61]の手法により生成する。この手法では、複雑な地表面・境界面形状を持った成層構造の入力データ（数値標高モデルによる深度データと物性値）に対し、バックグラウンド構造
格子を用いて、四面体要素とボクセル要素からなるハイブリッドメッシュを生成する。ボクセル要素は
更に 6 つの四面体要素に分割され、四面体要素は 10 節点の二次要素として用いる。[44][61]のメッシュ生成手法では、メッシュ生成対象領域の全域の層境界について、同じ水平方向のサイズのメッシュが生成される。局所的に更に高分解能なメッシュが必要な場合には、八分木分割パターンに従い、四面体、要素を一辺の長さが 1/2 の相似な四面体要素へと選択的に再分割することとする。
ここで、これらのメッシュ生成・有限要素計算プログラムにより計算される弾性応答が、メッシュサ
イズがある程度大きい場合、二次要素による離散化における数値誤差がメッシュサイズ ℎ に対する理論
的な収束速度である 𝑂(ℎ3 ) に従って収束することを、半無限均質媒質の弾性応答との比較によって確認する。グリーン関数ライブラリの計算設定を模した数値実験を設定し、これらの有限要素計算手法により計算される数値解の誤差の収束特性を、半無限弾性体における静弾性の変形応答の解析解[23]との比較により確認する。目的は、プログラムの veriﬁcation の一環ということに加え、メッシュサイズ・計算領
域サイズに対する数値誤差の定量化を行うための問題設定を決める助けとなる情報を得ることにある。
表 1.3.1 に問題設定と有限要素法による数値シミュレーションの設定を示す。日本列島地下の代表的な地震波速度パラメータと近い物性パラメータを持った半無限弾性体の深さ 20 km の部分に水平な平面断層を設定し、Section②-(b)で説明したものと空間的スケールの近い Bi-cubic B-spline 型のすべり分布（図 1.3.1）
を入力する。すべりをあたえる断層全体は 64km×64km の領域サイズとなる。16km×16km の格子サイズ
で重ね合わせた場合に Partition-of-unity 条件を満たすすべり分布となっており、最大すべり量は 0.444mである。解析解には岡田による変位量一様の矩形断層に対する応答の解析的表現 [23]を用いた。0.25- 34 -

km×0.25 km の矩形断層により B-spline 型すべり分布を離散化して計算したものを解析解ベースのReference solution とする。有限要素解析においては、計算領域中心部にすべりを入力し、計算領域端部の境界条件の数値解への影響を緩和するため、領域を広く設定する。そのうえで、𝑧 = 0 Kmの地表面における、領域中心部の 4km グリッドでの変位応答に注目する。計算用のメッシュについては、まず水平・鉛直方向ともに 8km 程度のサイズのメッシュを生成し、これをベースモデルとする。数値解を評価する
領域中心部において地表面から断層面までを含む、表 1.3.1 に示される mesh reﬁnement region を 1 段階、
2 段階、3 段階と再分割することで、メッシュサイズが 4km、2km、1km のメッシュを生成し、メッシュ
サイズに対する数値解の収束特性を評価する。変位応答を評価する全グリッドにおける変位三成分すべてを対象とした、参照解との誤差の L2 ノルムを評価する。図 1.3.2 に、メッシュサイズに対応する L2 ノルム誤差の収束の様子を示す。ここでは、地表面の−10km ≤ 𝓍 ≤ 10km 、−10km ≤ 𝑦 ≤ 10km の領域を評価対象としている。丸印で表されるメッシュサイズと誤差の関係が、点線で示される四面体二次要素における L2 ノルム誤差の理論的な収束速度である 𝑂(ℎ3 ) の直線上に載っていることがわかる。
ただし、ℎ はメッシュのサイズである。しかし、メッシュサイズが 1km となったところで、実際の誤差の減少量が理論的な収束速度から導かれるものより小さくなっている。これは、離散化による数値誤差が小さくなるにつれ、別の要因による数値誤差が顕在化して来ているものと思われる。別の主たる数値誤差要因
として計算領域端部の境界条件の影響を想定している。そこで、中心部が 1km メッシュのモデルの領域
を−880km ≤ 𝓍 ≤ 880km、−880km ≤ 𝑦 ≤ 880km、−980km ≤ 𝓏 ≤ 0km と拡張した設定で計算をしてみると、バツ印で表されるように数値誤差が拡張前と比べてかなり小さくなる。領域を更に広げていく
ことで、この数値誤差も更に小さくなり、理論値に近くなることが推測される。ただ、この時点で L2 ノルム誤差の変位応答 L2 ノルムに対する相対的な大きさは10−3 % 以下となっており、非常に小さい。以上のことから、対象とする実問題に近い空間スケールの問題において、モデル領域を十分広く確保した場
合に、メッシュサイズを 1km 程度とすることで十分数値誤差が小さくなることがわかる。この情報に基づき、Section③以降のグリーン関数ライブラリの数値誤差の議論において、最小メッシュサイズを 1kmとした計算モデルを基準とし、メッシュサイズとモデル領域が十分であるかを検討することが妥当だと判断する。またこの検討結果は、本プログラムによる数値誤差を論じる上で、メッシュサイズの影響と、変位固定の境界条件の影響の両方を議論する必要があることを支持する。

[50]は、グリーン関数ライブラリの生成の対象として、日本近海で巨大地震を起こすことが知られている南海トラフ域と日本海溝域を扱った。日本列島においては、長周期地震動を計算するための統一的な
地震波速度構造モデルとして JIVSM[24][25]が構築されている。この弾性構造データベースをもとに、南海トラフ域と日本海溝域それぞれについて有限要素モデルが生成され、計算に用いられた(図 1.3.3(a)(b))。
有限要素計算にはデカルト座標系を用いる一方、これらのデータは地理座標系で定義される。座標系の変換に際して、何らかの投影が行なわれることが多いが、モデル領域が大きくなるにつれ、投影による歪みの影響が無視できなくなる。これを避けるため、[50]では the International Terrestrial Reference Frameに従って the International Union of Geodesy and Geophysics と the Japanese Geodetic Datum において採用される Geodetic Reference System 1980 (GRS80)における座標系の定義と準拠楕円体に基づき、投影を伴わない変換が行なわれている。小断層としては、最大すべり量が 0.444m である 3 次の Bi-cubic B スプライン
形状(図 1.3.1)を持った単位すべりをプレート境界面上に配置した構造格子上の各点に入力する。すべり方- 35 -

向は以下の 2 方向を考慮する。1 方向は格子点におけるプレート境界面への接面に沈み込み方向を投影した方向とし、もう一方向は格子点におけるプレート境界面に対する法線ベクトルと沈み込みベクトルの双方に直行する方向とする。西南日本の南海トラフ域におけるプレート境界の主だった領域に Bi-cubic B
スプライン形状の小断層を導入するための格子点を配置すると、その個数は 460 個となり、沈み込み方
向とその垂直方向とを合わせて 920 個のグリーン関数を計算することとなった(図 1.3.3(c))。東北日本の日本海溝域に対しても同様に格子点を配置すると、その個数は 409 個となり、沈み込み方向とその垂直方向とを合わせて 818 個のグリーン関数を計算することとなった(図 1.3.3(d))。
メッシュ生成のために用いたパラメータは、日本海溝域での M9 地震に対する粘弾性応答を 1%以下の誤差レベル（応答レベルが小さな部分については期待される観測精度以下の誤差レベル）で計算するための設定を調べた[55]に基づいて設定した。これにより、層境界の水平方向要素サイズを 1km とした。
要素サイズ 1km は、グリーン関数ライブラリの計算設定と空間的スケールの近い問題での解析解との比較と収束速度の確認結果からも妥当と考えられる値である。この場合、鉛直方向のサイズはジオメトリに応じて変化する。また、[55]における最浅層の S 波速度が 3,550m/s であったのに比べ、今回使用した地震波速度構造においてはそれが大幅に小さい 700-2,000m/s 程度の層が存在したため、水平方向要素サイズが 1 km で十分であるかをメッシュサイズを変えた数ケースの事前解析により検証した。具体的には、上盤側において S 波速度が 700～2,000m/s 程度と小さい堆積層より浅い層における水平方向要素サイズを再分割により 0.5km、0.25km と変化させたメッシュを準備し、[55]で用いたものと同様の断層すべりを日本海溝域モデルに入力し、海域の観測点における弾性応答の変化が想定する観測誤差と比べて十分に小さい値として設定した 1％の範囲に収まる要素サイズを採用することとした。この結果、上盤側の堆積層
よ り 浅 い 層 の 水 平 方 向 メ ッ シ ュ サ イ ズ の み 0.5km と し た 。 メ ッ シ ュ 生 成 領 域 は 、2,500km×2,500km×1,100km と し て お り 、 そ の う ち 計 算 結 果 を 使 う 対 象 は 、 地 表 面 の 中 央 部
1,000km×1,000km の領域とした。この領域サイズは、計算結果を使う領域が[55]のケースより大幅に小さいことと、津波波源推定のための地殻変動計算における解析解との比較[62]等の過去の解析経験に基づき、[55]よりも水平方向に少し小さめに設定した。本研究の目的は、これらの計算モデルの設定が適切である
かどうかを、最新の計算機を用いてより網羅的に検証し、この設定に伴う数値誤差を定量化することである。
グリーン関数ライブラリにおいては、
地表面変位の値の多くが評価の対象となる。[50]においては、1,000個弱の入力断層すべり一つ一つに対し、モデル中央部の水平方向に 1,000km×1,000km の領域に対して、1km 間隔のグリッド(1,002,001 点)での変位が出力されている。そのすべてが公開対象となる。この全点での計算変位を、成分ごとに、ポイントワイズですべて評価することとする。Section②-(a)において有限要素法の数値解を解析解と比較したが、グリーン関数ライブラリで対象とする問題には解析解が当然存在しないため、同じアプローチは取れない。そこで、評価対象とするベースモデルを設定し、そこから見て相対的にメッシュが十分に細かい、あるいは領域が十分に広いモデルにおける数値解を参照解として用いる。このようなアプローチは、ベースモデルよりはるかに大きな計算を可能とする計算資源・計
算アルゴリズムを確保できて初めて可能となる。
参照解に対する計算値の評価指標として、絶対差と相対差の二つが考えられる。絶対差はεabsolute = |𝑢reference − 𝑢numerical|,
相対差は- 36 -

εrelative = |(𝑢reference − 𝑢numerical)/𝑢reference|のように定義する。ここで、ureference、unumerical はそれぞれ、参照解と評価対象の数値解を表す。絶対差は応答値の絶対値にも依存すると考えられるため、公平な比較のためには応答値によってスケーリングされた相対差を用いることが必要だと思われる。一方で、相対差は応答値が小さいときに値が大きくなる傾向にあるため、そのような場合は絶対差による評価が適当と考えられる。このように、両者は相補的な関係にあり、絶対差と相対差の両方を評価する必要がある。また 1,000 ケース弱における全 1,002,001
点に対するポイントワイズの評価は、計算という観点からは困難でないが、結果を可視化して確かめるという意味では全く容易ではない。差の情報を縮約して確認するための評価方法が必要となる。このことを踏まえ、本研究では絶対差と相対差の閾値を設定した上で、以下のような方法で評価を行う。まず
地表面の全 1,002,001 点の各点に対して相対差を計算し、相対差の閾値に収まるかどうかを評価する。相対差が閾値を超える点について絶対差を評価する。絶対値の閾値を超える点の全 1,002,001 点に対する割合を評価する。この割合を用いた指標 EI1 を数式で表現すると
EI1 = card(𝐴)/𝑁,
ただし
thres
thres
∶= {𝑃𝑖 |1 ≤ 𝑖 ≤ 𝑁, 𝜀relative(𝐱𝑃𝑖 ) > 𝜀relative∧ 𝜀absolute(𝐱𝑃𝑖 ) > 𝜀absolute}𝐴となる。ここで 𝑁、 𝑖、 𝑃𝑖 、xPi はそれぞれ、評価点の数（1,002,001）、地表面上の評価点のインデック
ス、i 番目の評価点、𝑃𝑖 の座標である。このような指標を導入することで、参照解とベースモデルによる数値解の差を、絶対差と相対差の両方を考慮した形で集約して評価することができる。また、相対差の閾値を超えたが絶対値の閾値を超えなかった点での絶対値の最大値も評価する。この最大値を用いた指標 EI2 を数式で表現するとEI2 = max 𝐵,
ただし𝐵∶= {𝜀absolute(𝐱𝑃𝑖 ) | 𝑃𝑖 ∈ 𝐶},

𝐶

thres
thres
∶= {𝑃𝑖 |1 ≤ 𝑖 ≤ 𝑁, 𝜀relative(𝐱𝑃𝑖 ) > 𝜀relative∧ 𝜀absolute(𝐱𝑃𝑖 ) < 𝜀absolute}となる。この量を調べることにより、応答の差が閾値を超えることがなかった小断層についても、性質をある程度把握することができる。これらの評価は各小断層に対する応答について実施し、変位の各 3
thres
thres
成分について独立に行う。 𝜀absoluteと𝜀relative
は、目的に応じた十分小さい値に設定する必要がある。本
thres
thres
研究では、𝜀absolute= 10−4 (m)と𝜀relative= 10−3 を基本の設定とする。これらの閾値を客観的に決めるのは難しいが、過去の観測記録と観測の誤差レベルと照らし合わせることで、現行の地震学において期待される用途に対して十分な閾値であると判断した。
本研究では、[50]で発表したグリーン関数ライブラリを評価することを目指しているので、ベースモデ
ルはこの研究で用いたものと基本的に同様とする。しかし、そこで使用された日本列島の有限要素モデ- 37 -

ルは、水平方向の最小メッシュサイズを 1km としてメッシュ生成したのち、地表面付近の四面体要素を
さらに 0.5km に再分割したものであった。ここでは、要素再分割の影響を評価しなおしたいので、再分
割を施さずに最小メッシュサイズを 1km としたままのモデルをベースモデルとする。また今回の計算では、地表面とプレート境界面に位置する節点を、それらの面の形状が入力の DEM データにより忠実になるように、メッシュ生成の後に移動する処理を行っている。
評価対象とするベースモデルに対して、要素再分割によりメッシュを細かくした計算モデルを作成する。ここでは、四面体要素を 2 段階分割し、最小要素サイズを 0.25km としたモデルを考える。このモデ
ルを Reﬁned2 モデルと呼ぶ。
2 段階分割によりベースモデルに対するメッシュサイズは(1/2)2 となるため、二次要素のメッシュサイズに対する理論的な収束速度𝑂(ℎ3) を考慮すると、離散化による数値誤差は(1/2)6 となる。つまり 2 段階再分割モデルとベースモデルの差の評価により、数値誤差の 1−(1/2)6 = 63/64を評価できていることになるため、両者の差を数値誤差の近似値として十分使用できると考えられる。
そこで、𝑢reference = 𝑢Reﬁned2 と考えて評価をおこなう。ここで、𝑢Reﬁned2 は Reﬁned2 モデルを用いて得られた数値解を表す。また、メッシュサイズの変化に対する数値誤差のふるまいを確認するため、1 段階分
割による最小メッシュサイズ 0.5km のモデル、Reﬁned1 モデルに関しても、参照解の差を併せて評価する。
要素の再分割は、変位を評価するモデル中央部の地表面付近と断層すべりを入力するプレート境界を含む領域に対して行う(図 1.3.4)。
注意するべきことは、この要素再分割においては、計算モデルのジオメトリは固定としていることである。再分割により新しくできた節点位置は、分割前の要素の辺上に存在させ、そこからの移動は行わない。このようにする理由は、形状情報として使っている JIVSM の DEM データに、1km 以下のスケー
ルの詳細な情報が多く含まれていないからである。もっとも、Section ②-(a)にて説明したメッシュ生成
アルゴリズムでは、選択した最小メッシュサイズに依存した形状の加工が行われる。計算モデルを作成する際の入力データからの形状の変更は、この手法に限らず数値シミュレーション一般に生じうる過程である。この形状加工の影響が数値解に及ぼしている影響は、今回のアプローチでは評価できない。このような場合に手続きの透明性を担保するには、地表面とプレート境界の形状が入力の DEM データにより忠実になるようにして形状加工を小さくとどめるとともに、形状加工後の形状もデータとしてアクセス可能とすることが考えられる。同時に、形状加工の影響を評価する手法も、今後の課題として検討の必要がある。
評価対象とするベースモデルに対して、計算領域を広げた計算モデルを作成する。ここでは、x 方向、
y 方向、z 方向すべての方向に 500km 広げたモデル、Wide モデルを用いることとする(図 1.3.4)。つまり、
𝑢reference = 𝑢Wide と考えて評価をおこなう。ここで、𝑢Wide は Wide モデルを用いて得られた数値解を表す。計算領域サイズに関する参照モデルの作成に際しても、メッシュサイズに関する検討と同様に、理論的な数値誤差の予測に基づいて決められることが望ましい。しかし、計算領域端部の境界条件設定の及ぼす数値誤差と計算領域サイズの関係性はケースバイケースと考えられ、理論的関係を考えることが難しい。一方、のちの検討において示されるように、計算領域サイズによる数値誤差の影響は、メッシ
ュサイズによる影響に比べて相対的に小さいため、計算領域サイズに関する検討についてはこのような簡便な比較で十分だとみなすこととする。
図 1.3.5 と図 1.3.6 に、Base モデルと Reﬁned2 モデルの比較における、各小断層におけるEI1 とEI2の値- 38 -

のマップを変位成分ごとに示す。ここでは、沈み込みと反対方向のすべりに対する応答のみ示した。南海トラフモデルの伊豆半島北部の地下に置いたものを除くすべての小断層に対し、地表面の変位評価点
1,002,001 点のうち数値誤差の閾値を超える点の数は多くとも 0.01%以下、つまり 100 個以下のオーダーに抑えられている。閾値を超える点を持つ小断層はトラフ・海溝軸近傍の小断層に限られており、それらからある程度離れた小断層においては、すべての変位評価点において数値誤差が閾値を下回っている。
閾値を下回った点の中で最大の絶対差が大きくなっているのもトラフ・海溝軸近傍の小断層であり、ある程度離れた小断層では、最大の絶対差は 0.1mm 以下のオーダーに抑えられている。これより、トラフ・海溝軸近傍の小断層に対する応答の使い方に注意をすれば、ベースモデルの離散化に関する数値誤差は要求する水準を十分に満たしていることが分かった。なお、Reﬁne1 モデルと Reﬁned2 モデルについても同様の比較を行なったところ、Base モデルと Reﬁned2 モデルの比較に比べて、差がさらに小さくなっていることが分かった（図は省略）。これらの比較においてはメッシュサイズが十分細かくなっていると考えられるため、図 1.3.2 に示した数値実験の結果と同様に、収束速度に関する議論は当てはまらないと考えられる。だが、メッシュを細かくすることで数値解が Reﬁne2 により近づいていることがわかり、数値解が期待通りのふるまいをしていることが確認できた。
図 1.3.7 と図 1.3.8 に、Base モデルと wide モデルの比較における、各小断層におけるEI1 とEI2の値のマップを x 変位成分について示す。ここでは、沈み込みと反対方向のすべりに対する応答のみ示した。南海トラフモデル、日本海溝モデルともに、地表面の変位評価点 1,002,001 点のうち数値誤差の閾値を超える点は存在せず、最大の絶対差もカラーバーでは判別できないほど小さいレベルに抑えられている。これより、ベースモデルの計算領域サイズに関する数値誤差は要求する水準を十分に満たしており、離散化に関する数値誤差と比べてもマイナーな影響しかないことが分かった。
表 1.3.2 に、数値誤差の定量化のために用いた各計算モデルの自由度を示す。南海トラフモデル・日本
海溝モデルともに、Reﬁned2 モデルは 400 億自由度程度となっており、ベースモデルの 8-10 倍程度の規模となっている。必要な計算資源も当然 Reﬁned2 モデルが突出している。例えば、日本海溝については、
IO 等計算以外にかかった時間も含めた参考情報であるが、Base モデルは「富岳」の 512 ノードを使って8.5 時間で終了したのに対し、Reﬁned2 モデルに対しては「富岳」の 4,096 ノードを使って 31.5 時間かかっている。[50]は、グリーン関数ライブラリの計算に JCAHPC[1]の Oakforest-PACS[56]（1 ノード辺り理
論性能 3TFLOPS 程度）を用いた。必要だったのは 512 ノードを用いて 16 時間余りであり、今回の Base
モデルもそれと同様のオーダーとなる。一方で、Reﬁned2 モデルに必要な計算は、Oakforest-PACS の約半分を同じ時間占有することに相当する。これは本機の通常の運用形態においては実行不可能な計算規模である。
「富岳」のような巨大なマシンが使えるようになる恩恵は、全系を使ったフラッグシップ的な計算だけではなく、多数のグリーン関数の数値誤差定量化のようなより実用的かつ多少の試行錯誤が必要となる検討に対しても波及すると言える。
③の要点をまとめる。地震学・ハザード評価に広い用途のあるプレート境界断層すべりに対する地表
面変位のグリーン関数ライブラリに対し、数値誤差の定量化を行う枠組みを提案した。この枠組みは、
よりメッシュを細かくした・領域を広くしたモデルを参照解として設定し、そこからのずれを数値誤差として扱うものである。その数値誤差を、相対差と絶対差の両面から評価し、設定した閾値との関係を 2
つの評価指標にまとめ、対象とする地域・変位成分・入力すべり成分ごとに可視化を行う。この枠組み- 39 -

を、[50]の南海トラフ域・日本海溝域のグリーン関数ライブラリへと適用した。よりメッシュを細かくした・領域を広くしたより大規模な計算モデルでの計算実行は、高速かつスケーラブルな弾性応答計算ア
ルゴリズムを活用して「富岳」を用いることにより可能となった。可視化の結果、誤差が閾値を超えた
地点は、食い違い変位モデルが入力されている地表に露出した断層面の部分に相当するトラフ軸と海溝軸の近傍にほぼ限られていることがわかった。メッシュサイズの異なる 2 つのモデルで計算したグリーン関数を用いて、史上最大の振幅変位データを持つ 2011 年東北地方太平洋沖地震の断層すべり推定問題の解を比較した。2 つのすべり推定結果の間に非常に良い一致が見られたことから、提案した誤差指標と閾値で検証したグリーン関数ライブラリは、現在地震学で期待されている多くのアプリケーションに対して十分な性能・精度を持っていると考えられる。一方で、トラフ軸や海溝軸のごく近傍で観測が行われる状況下では、提案した誤差の閾値やグリーン関数ライブラリの計算設定の修正が必要となる場合がある。本研究で提案した数値誤差評価手法は、他の沈み込み帯や断層を対象としたグリーン関数ライブ
ラリにも適用できることはもちろんのこと、断層面上のせん断応力計算等、他の物理量の数値誤差の評価にも適用可能である。また、グリーン関数ライブラリの形であるかにかかわらず、弾性応答の計算に対して誤差の影響がクリティカルだと推測される部分の数値誤差の定量化を行う際に、選択的に提案のアプローチを適用することも有用と考える。
表 1.3.1 収束特性確認のための計算設定。𝒄(𝒙, 𝒚) は図 10 に示す双 3 次 B スプラインの関数形状に従った空間すべり分布を特徴づける関数である。
Half-space

Fault plane

Dislocation

FEM model region

Mesh reﬁnement region

Surface: 𝓏 = 0 km
Rigidity: µ = 31.2 GPa
Poisson’s ratio: ν = 0.238-32 km ≤ 𝓍 ≤ 32 km-32 km ≤ 𝑦 ≤ 32 km
𝓏 = -20 km∆𝑑𝑥 = 0.819 × 𝑐(𝑥, 𝑦) m
∆𝑑𝑦 = −0.534 × 𝑐(𝑥, 𝑦) m
∆𝑑𝑧 = 0 m
-640 km ≤ x ≤ 640 km-640 km ≤ 𝑦 ≤ 640 km-500 km ≤ 𝓏 ≤ 0 km
-160 km ≤ x ≤ 160 km-160 km ≤ 𝑦 ≤ 160 km-30 km ≤ 𝓏 ≤ 0 km

- 40 -

図 1.3.1 𝒄(𝒙, 𝒚)の空間分布

図 1.3.2 メッシュサイズに対する誤差の収束の様子。丸印が標準の計算領域、バツ印が拡張した計算領域によるモデルの結果のプロット。点線は理論的な収束速度である𝑶(𝒉𝟑 )を表す。

- 41 -

図 1.3.3 (a)(b) 南海トラフと日本海溝の沈み込み帯の FE モデルの概要。(c)(d) 南海トラフと日本
海溝の沈み込み帯の GF 計算対象領域の概観。灰色の点は、サブ断層に適用した B-スプライン関数型のすべり分布の中心点を示す。三角形は GEONET、逆三角形は測地観測の代表的な海底観測
点をそれぞれ示している。図は[50]に掲載のものを改変して再掲したものである。

- 42 -

図 1.3.4 数値誤差の推定を目的として、(a)南海トラフ、(b)日本海溝の FE モデルに施した修正の模式図。白色文字・記号は Reﬁned1・Reﬁned2 でのメッシュ再分割領域、灰色文字・記号は、Wide
モデルでのモデル領域の拡大をそれぞれ示す。

図
図 1.3.5
1.3.6 南海トラフ海域において、
日本海溝域において、沈み込み方向と逆の方向に断層がすべる場合の変位𝔁,
沈み込み方向と逆の方向に断層が滑った場合の変位𝔁,𝒚, 𝒚,𝔃の各成分につい𝔃の各成分につ- 43 ) を表している．白色
いての、Baseての、Base モデルとモデルとReﬁned2Reﬁned2モデルの比較。左図の円の色は各小断層の𝐥𝐨𝐠モデルの比較。左図の円の色は各小断層の𝐥𝐨𝐠
𝟏𝟎 (𝐄𝐈
𝟏𝟎 (𝐄𝐈
𝟏 ) 𝟏を表している．白色は
は𝐄𝐈
𝐄𝐈𝟏 =
𝟎 𝟎、黒色は𝐄𝐈、黒色は𝐄𝐈
𝟎. 𝟏%を示す．右の円の色は、各小断層に対する𝐄𝐈
𝟎. 𝟏%を示す．右の円の色は、各小断層に対する𝐄𝐈
ID
𝟏 =
𝟏 >
𝟏 >
𝟐 を表す。円内の数字は小断層𝟐 を表す。円内の数字は小断層
ID
を表す。
を表す。

図 1.3.7 南海トラフ海域において、沈み込み方向と逆の方向に断層が滑った場合の変位𝔁, 𝒚, 𝔃の各成分についての、Base モデルと Wide モデルの比較。左図の円の色は各小断層の𝐥𝐨𝐠𝟏𝟎 (𝐄𝐈𝟏 ) を表している．白色は𝐄𝐈𝟏 = 𝟎 、黒色は𝐄𝐈𝟏 > 𝟎. 𝟏%を示す．右の円の色は、各小断層に対する𝐄𝐈𝟐 を表す。円内の数字は小断層
ID を表す。

図 1.3.8 日本海溝域において、沈み込み方向と逆の方向に断層がすべる場合の変位𝔁, 𝒚, 𝔃の各成分についての、Base モデルと Wide モデルの比較。左図の円の色は各小断層の𝐥𝐨𝐠𝟏𝟎 (𝐄𝐈𝟏 ) を表している．白色は𝐄𝐈𝟏 = 𝟎 、黒色は𝐄𝐈𝟏 > 𝟎. 𝟏%を示す．右の円の色は、各小断層に対する𝐄𝐈𝟐 を表す。円内の数字は小断
層 ID を表す。

- 44 -

表 1.3.2 数値誤差の定量化に使用した各計算モデルの自由度
Model
Base
Reﬁned1
Reﬁned2
Wide

DOF (billion)
Nankai trough Japan trench
5.2
4.0
9.0
8.7
39.4
39.9
7.4
7.0

[1] Hisada, Y., and S. Yamamoto. "One-, two-, and three-dimensional site effects in sediment-filled basins."
Proceedings of the 11th world conference on earthquake engineering. 1996.
[2] T . Ichimura, M. Hori, J. Bielak A hybrid multiresolution meshing technique for finite element
three-dimensional earthquake ground motion modelling in basins including topography, Geophysical JournalInternational, Volume 177, Issue 3, Pages 1221 1232, June 2009.
[3] Ichimura, T., Fujita, K., Tanaka, S., Hori, M., Lalith, M., Shizawa, Y., & Kobayashi, H. (2014, November).
Physics-based urban earthquake simulation enhanced by 10.7 BlnDOF× 30 K time-step unstructured FE
non-linear seismic wave simulation. In SC'14: Proceedings of the International Conference for High
Performance Computing, Networking, Storage and Analysis (pp. 15-26). IEEE.
[4] Ichimura, T., Fujita, K., Quinay, P. E. B., Maddegedara, L., Hori, M., Tanaka, S., Shizawa, Y., Kobayashi, H.
& Minami, K. (2015, November). Implicit nonlinear wave simulation with 1.08 T DOF and 0.270 T
unstructured finite elements to enhance comprehensive earthquake simulation. In SC'15: Proceedings of theInternational Conference for High Performance Computing, Networking, Storage and Analysis (pp. 1-12).
IEEE.
[5] Bielak, J., Ghattas, O., & Kim, E. J. (2005). Parallel octree-based finite element method for large-scale
earthquake ground motion simulation. CMES - Computer Modeling in Engineering and Sciences, 10(2), 99–
112. https://doi.org/10.3970/cmes.2005.010.099[6] Melosh, H., & Raefsky, A. (1981). A simple and efficient method for introducing faults into finite element
computations. Bulletin of the Seismological Society of America, 71(5), 1391–1400.
https://doi.org/10.1785/BSSA0710051391[7] Lysmer, J., & Kuhlemeyer, R. L. (1969). Finite dynamic model for infinite media. Journal of the Engineering
Mechanics Division, 95(4), 859-878.
[8] Cerjan, C., Kosloff, D., Kosloff, R., & Reshef, M. (1985). A nonreflecting boundary condition for discrete
acoustic and elastic wave equations. Geophysics, 50(4), 705–708. https://doi.org/10.1190/1.1441945[9] Hisada, Yoshiaki. "An efficient method for computing Green's functions for a layered half-space with sources
and receivers at close depths." Bulletin of the Seismological Society of America 84.5 (1994): 1456-1472.
[10] Anderson, J. G. (2004). Quantitative measure of the goodness-of-fit of synthetic seismograms. Proceedings ofthe 13th World Conference on Earthquake Engineering, 243, 243.
http://www.iitk.ac.in/nicee/wcee/article/13_243.pdf- 45 -

[11] Wu, Z. L., & Chen, Y. T. (2003). Definition of seismic moment at a discontinuity interface. Bulletin of the
Seismological Society of America, 93(4), 1832-1834.
[12] 林田智宏, 岡元太郎, 竹中博士: 速度-応力型スタガード格子差分法を用いた 2 次元及び 3 次元地震波
動計算コードの作成, 九州大学理学部研究報告 地球惑星科学, 20 巻, 3 号, pp. 99-110, 1999.
[13] Graves, R.: Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences,
Bull. Seism. Soc. Am., Vol. 86, No. 4, pp. 1091-1106, 1996.
[14] Clayton, R., Engquist, B.: Absorbing boundary conditions for acoustic and elastic wave equations, Bull. Seism.
Soc. Am., Vol. 67, No. 6, pp. 1529-1540, 1977.
[15] Zeng, C., Xia, J., Miller, R., Tsoflias, G.: An improved vacuum formulation for 2D finite-difference modelingof Rayleigh waves including surface topography and internal discontinuities, Geophysics, Vol. 77, No. 1, pp.
T1-T9, 2012.
[16] Pitarka A.: 3D Elastic finite-difference modeling of seismic motion using staggered grids with nonuniform
spacing, Bull. Seismol. Soc. Am., Vol. 89, No. 1, pp. 54-68, 1999.
[17] 笠松健太郎, 山中浩明: 中小地震の地震動シミュレーションによる 2004 年新潟県中越地震震源域の深部地盤構造モデルの構築, 物理探査, 59 巻, 5 号, pp. 475-484, 2006.
[18] 笠松健太郎, 加藤研一: 大規模並列計算による 2011 年東北地方太平洋沖地震（Mw9.0）の長周期地震
動シミュレーション, 地球惑星科学連合大会, 2016.
[19] 笠松健太郎, 渡辺哲史, 野澤貴, 加藤研一, 酒井慎一: 地形と海水の影響を考慮した 2009 年伊豆半島沖の地震 M5.1 の地震動シミュレーション, 日本建築学会大会学術講演梗概集, pp. 149-150, 2020.
[20] 地震調査研究推進本部: 「長周期地震動予測地図」２０１２年試作版, 2012.
https://www.jishin.go.jp/evaluation/seismic_hazard_map/lpshm/12_choshuki/[21] 上野寛, 畠山信一, 明田川保, 舟崎淳, 浜田信生 : 気象庁の震源決定方法の改善－浅部速度構造と重み関数の改良－, 験震時報, 65, pp. 1507-1520, 2002.
[22] Y Okada. Surface deformation due to shear and tensile faults in a half-space. Bulletin of the SeismologicalSociety of America, 75(4):1135–1154, 1985.
[23] Yoshimitsu Okada. Internal deformation due to shear and tensile faults in a half-space. Bulletin of the
Seismological Society of America, 82(2):1018–1040, 1992.
[24] Kazuki Koketsu, Hiroe Miyake, Yasuhisa Tanaka, et al. A proposal for a standard procedure of modeling 3-D
velocity structures and its application to the Tokyo metropolitan area, Japan. Tectonophysics, 472(1- 4):290–
300, 2009.
[25] Kazuki Koketsu, Hiroe Miyake, and Haruhiko Suzuki. Japan integrated velocity structure model version 1.
Proceedings of the 15th World Conference on Earthquake Engineering, (1773), 2012.
[26] Gabi Laske, Guy Masters, Zhitu Ma, and Mike Pasyanos. Update on CRUST1.0 – A 1–degree global model ofEarth’s crust. EGU General Assembly Vienna, Austria, 15:2658, 2013.
[27] Ayako Nakanishi, Narumi Takahashi, Yojiro Yamamoto, Tsutomu Takahashi, Seckin Ozgur Citak, Takeshi
Nakamura, Koichiro Obana, Shuichi Kodaira, and Yoshiyuki Kaneda. Three-dimensional plate geometry andP-wave velocity models of the subduction zone in SW Japan: Implications for seismogenesis. Geology and
Tectonics of Subduction Zones: A Tribute to Gaku Kimura, 534:69, 2018.

- 46 -

[28] Shinichi Miyazaki and Hatanaka Yuki. The outlines of the GEONET (in Japanese). Meteorol. Res. Note,
192:105–131, 1998.
[29] C. David Chadwell and Fred N. Spiess. Plate motion at the ridge-transform boundary of the south Cleftsegment of the Juan de Fuca Ridge from GPS-Acoustic data. Journal of Geophysical Research,
113(B4):B04415, apr 2008.
[30] Mariko Sato, Tadashi Ishikawa, Naoto Ujihara, Shigeru Yoshida, Masayuki Fujita, Masashi Mochizuki, andAkira Asada. Displacement above the hypocenter of the 2011 Tohoku-Oki earthquake. Science (New York,N.Y.), 332(6036):1395, 2011.
[31] T Iinuma, R Hino, M Kido, D Inazu, Y Osada, Y Ito, M Ohzono, H Tsushima, S Suzuki, H Fujimoto, and S
Miura. Coseismic slip distribution of the 2011 oﬀ the Paciﬁc Coast of Tohoku Earthquake (M9.0) reﬁned bymeans of seaﬂoor geodetic data. Journal of Geophysical Research, 117(B7):1–18, Jul 2012.
[32] Fumiaki Tomita, Motoyuki Kido, Yukihito Osada, Ryota Hino, Yusaku Ohta, and Takeshi Iinuma. First
measurement of the displacement rate of the paciﬁc plate near the Japan trench after the 2011 Tohoku-oki
earthquake using gps/acoustic technique. Geophysical Research Letters, 42(20):8391–8397, 2015.
[33] Yusuke Yokota, Tadashi Ishikawa, Shun-ichi Watanabe, Toshiharu Tashiro, and Akira Asada. Seaﬂoor geodetic
constraints on interplate coupling of the Nankai Trough megathrust zone. Nature, 534(7607):374, 2016.
[34] Y Kaneda, K Kawaguchi, E Araki, A Sakuma, H Matsumoto, T Nakamura, S Kamiya, K Ariyoshi, T Baba, M
Ohori, et al. Dense Ocean ﬂoor Network for Earthquakes and Tsunamis (DONET)- Development and Data
application for the mega thrust earthquakes around the Nankai trough. 1:1453, 2009.
[35] Yoshihiro Ito, Takeshi Tsuji, Yukihito Osada, Motoyuki Kido, Daisuke Inazu, Yutaka Hayashi, Hiroaki
Tsushima, Ryota Hino, and Hiromi Fujimoto. Frontal wedge deformation near the source region of the 2011
Tohoku-oki earthquake. Geophysical Research Letters, 38(7), 2011.
[36] M. Kido, Y. Osada, H. Fujimoto, R. Hino, and Y. Ito. Trench-normal variation in observed seaﬂoor displacements associated with the 2011 Tohoku-Oki earthquake. Geophysical Research Letters, 38, L24303, Dec
2011.
[37] Y Kaneda, K Kawaguchi, E Araki, H Matsumoto, T Nakamura, S Kamiya, K Ariyoshi, T Hori, T Baba, and N
Takahashi. Development and application of an advanced ocean ﬂoor network system for megathrustearthquakes and tsunamis. In Seaﬂoor observatories, pages 643–662. Springer, 2015.
[38] Didier Massonnet, Marc Rossi, Ce´sar Carmona, Fre´de´ric Adragna, Gilles Peltzer, Kurt Feigl, and Thierry
Rabaute. The displacement ﬁeld of the landers earthquake mapped by radar interferometry. Nature,364(6433):138, 1993.
[39] TJ Wright, BE Parsons, JA Jackson, M Haynes, EJ Fielding, PC England, and PJ Clarke. Source parameters of
the 1 October 1995 dinar (turkey) earthquake from sar interferometry and seismic bodywave modelling. Earthand Planetary Science Letters, 172(1-2):23–37, 1999.
[40] Tomokazu Kobayashi. Earthquake rupture properties of the 2016 Kumamoto earthquake foreshocks (M j 6.5
and M j 6.4) revealed by conventional and multiple-aperture InSAR. Earth, Planets and Space, 69(1):7, 2017.
[41] Timothy Masterlark. Finite element model predictions of static deformation from dislocation sources in a
subduction zone: Sensitivities to homogeneous, isotropic, Poisson-solid, and half-space assumptions. Journal- 47 -

of Geophysical Research, 108(B11):2540, 2003.
[42] T. Ichimura, R. Agata, T. Hori, K. Hirahara, and M. Hori. Fast numerical simulation of crustal deformation
using a three-dimensional high-ﬁdelity model. Geophysical Journal International, 195(3):1730–1744, Sep
2013.
[43] C. Kyriakopoulos, T. Masterlark, S. Stramondo, M. Chini, and C. Bignami. Coseismic slip distribution for the
Mw 9 2011 Tohoku-Oki earthquake derived from 3-D FE modeling. Journal of Geophysical Research: Solid
Earth, 118(February):3837–3847, Jul 2013.
[44] Tsuyoshi Ichimura, Ryoichiro Agata, Takane Hori, Kazuro Hirahara, Chihiro Hashimoto, Muneo Hori, andYukitoshi Fukahata. An elastic/viscoelastic ﬁnite element analysis method for crustal deformation using a 3-D
island-scale high-ﬁdelity model. Geophysical Journal International, 206(1):114–129, 2016.
[45] Charles A Williams and Laura M Wallace. Eﬀects of material property variations on slip estimates for
subduction interface slow-slip events. Geophysical Research Letters, 42(4):1113–1121, 2015.
[46] Hom Nath Gharti, Leah Langer, and Jeroen Tromp. Spectral-inﬁnite-element simulations of coseismic andpost-earthquake deformation. Geophysical Journal International, 216(2):1364–1393, 2018.
[47] Charles A Williams and Laura M Wallace. The impact of realistic elastic properties on inversions of shallow
subduction interface slow slip events using seaﬂoor geodetic data. Geophysical Research Letters,
45(15):7462–7470, 2018.
[48] Leah Langer, Hom Nath Gharti, and Jeroen Tromp. Impact of topography and three-dimensional heterogeneity
on coseismic deformation. Geophysical Journal International, 217(2):866–878, 2019.
[49] Huihong Cheng, Bei Zhang, Luyuan Huang, Huai Zhang, and Yaolin Shi. Calculating coseismic deformationand stress changes in a heterogeneous ellipsoid earth model. Geophysical Journal International, 216(2):851–
858, 2019.
[50] Takane Hori, Ryoichiro Agata, Tsuyoshi Ichimura, Kohei Fujita, Takuma Yamaguchi, and Takeshi Iinuma.
High-ﬁdelity elastic Green’s functions for subduction zone models consistent with the global standard geodetic
reference system. Earth, Planets and Space, 73(1):1–12, 2021.
[51] The American Society of Mechanical Engineers. ASME V&V 10-2006 - Guide for Veriﬁcation and Validationin Computational Solid Mechanics. The American Society of Mechanical Engineers, 2006.
[52] Jacobo Bielak, Robert W Graves, Kim B Olsen, Ricardo Taborda, Leonardo Ram´ırez-Guzma´n, Steven M
Day, Geoﬀrey P Ely, Daniel Roten, Thomas H Jordan, Philip J Maechling, et al. The shakeout earthquake
scenario: Veriﬁcation of three simulation sets. Geophysical Journal International, 180(1):375–404, 2010.
[53] Ruth A Harris, Michael Barall, R Archuleta, E Dunham, B Aagaard, Jean Paul Ampuero, Harsha Bhat, V
Cruz-Atienza, L Dalguer, Phillip Dawson, et al. The SCEC/USGS dynamic earthquake rupture codeveriﬁcation exercise. Seismological Research Letters, 80(1):119–126, 2009.
[54] Brittany A Erickson, Junle Jiang, Michael Barall, Nadia Lapusta, Eric M Dunham, Ruth Harris, Lauren S
Abrahams, Kali L Allison, Jean-Paul Ampuero, Sylvain Barbot, et al. The community code veriﬁcation
exercise for simulating sequences of earthquakes and aseismic slip (SEAS). Seismological Research Letters,
91(2A):874–890, 2020.
[55] Ryoichiro Agata, Tsuyoshi Ichimura, Kazuro Hirahara, Mamoru Hyodo, Takane Hori, Chihiro Hashimoto, and- 48 -

Muneo Hori. Numerical Veriﬁcation Criteria for Coseismic and Postseismic Crustal Deformation Analysis
with Large-scale High-ﬁdelity Model. Procedia Computer Science, 51:1534–1544, 2015.
[56] [Online]. JCAHPC, https://www.cc.u-tokyo.ac.jp/supercomputer/ofp/service/[57] [Online]. Oakforest-PACS, http://www.cc.u-tokyo.ac.jp/system/ofp/index-e.html.
[58] William L Oberkampf and Timothy G Trucano. Veriﬁcation and validation in computational ﬂuid dynamics.
Progress in aerospace sciences, 38(3):209–272, 2002.
[59] HJ Melosh and A Raefsky. A simple and eﬃcient method for introducing faults into ﬁnite element
computations. Bulletin of the Seismological Society of America, 71(5):1391–1400, 1981.
[60] Kohei Fujita, Tsuyoshi Ichimura, Kentaro Koyama, Hikaru Inoue, Muneo Hori, and Lalith Maddegedara. Fastand scalable low-order implicit unstructured ﬁnite-element solver for earth’s crust deformation problem. In
Proceedings of the Platform for Advanced Scientiﬁc Computing Conference, article No.11, pages 1–10. ACM,
2017.
[61] Tsuyoshi Ichimura, Muneo Hori, and Jacobo Bielak. A hybrid multiresolution meshing technique for ﬁnite
element three-dimensional earthquake ground motion modelling in basins including topography. GeophysicalJournal International, 177(3):1221–1232, 2009.
[62] [Online].「日本海溝・千島海溝沿いの 巨大地震モデルの検討について」
https://www.bousai.go.jp/jishin/nihonkaiko_chishima/model/pdf/sankozuhyou.pdf（２）地震に関する災害被害予測のための大規模アプリケーションの開発
本項目では、計算科学・計算機科学の最先端技術を創出しこれを駆使することで、地震に関する災害
被害予測のためのアプリケーション群の基幹数値解析コンポーネントである、地震・地震動・地殻変動計算等を行う大規模数値解析用コード群に対し、「富岳」実機上でより良好な性能を発揮するように研究開発を行っている。
本年度は、「富岳」の計算機構の特性を考慮した上で、「富岳」の性能を引き出すように、solver の前処理等の改良を行い、アプリケーション群の性能向上を図った。また、超大規模解析モデルを構築し、これを用いた解析を行うことで、解析手法の有効性を示した。具体的には、「①シミュレーション内で生成されるデータの学習によるシミュレーション性能の向上及び超高詳細断層-構造物系モデルを用いた地震シミュレーション」、「②シミュレーション内で生成されるデータの学習による反復 solver 性能の向上」、「③ニューラルネットワークによりデータ構造を組み換えた効率的な前処理の高度化」、「④開発された大規模数値解析コード群の実展開を踏まえた解析手法の安定化及び有効性の確認」を行った。なお、①は The International Conference on High Performance Computing in Asia-Pacific Region (HPC Asia 2022)にてBest paper award を、③は Eighth Workshop on Accelerator Programming Using Directives (WACCPD) @SC21
にて Honorable mention を、④は HPC Asia 2022 Best Student Poster Award in Memory of Hiroshi Nakashima の
1st Prize をいただいた。以下に、上記の①～④の具体的な成果をまとめる。
①シミュレーション内で生成されるデータの学習によるシミュレーション性能の向上及び超高詳細断層構造物系モデルを用いた地震シミュレーション- 49 -

①-(a) はじめに都市・社会への災害の脅威が高まるなか、post/with コロナ時代へ向けて digital transformation を踏まえた都市・社会の急激な変容がさらに加速されている。過去から大きく変容する状況において、過去の被害経験に頼るだけでは十分な被害軽減を期待することは難しい。
一方で、都市・社会の digital transformation
により、我々はいままでよりもさらに高詳細な都市・社会データを手に入れつつある。過去の被害経験と高詳細都市・社会データによるシミュレーションを組み合わせることで、いままでとは一線を画す信頼性の被害軽減策実現が期待されるが、この解析は巨大かつ人類が未経験のものであり、この実現にはSupercomputing の貢献が必須となる。
地震は、i)断層破壊により生成された地震波が地下深部を伝播し、ii)地表近傍の堆積層でその波動が増幅され、iii)増幅された地震波により都市・社会・構造物が強く揺らされるという非線形動的問題である（図 2.1.1）。従前から、このような社会的要請が非常に強い地震問題への Supercomputing の貢献が強く期待されてきており、様々な成果が創出されている(e.g. [1])。例えば、i)のための手法として[2][3][4][5]が、ii)のための手法として[6][7]が、ii)+iii)のための手法として[8]などが Supercomputing Conference で発表されており、地震のサイエンス及び震災軽減のための解析能力向上への大きな貢献が期待されている。

図 2.1.1 地震の諸過程。i) 断層破壊により生じた波の地下深部での波動伝播、ii) 地表近傍の堆積層での波動増幅、iii) 都市・社会・構造物が強く揺らされる。本来であれば、これら i)ii)iii)の一括解析が望まれるが、膨大な解析コストからその実現は難しいとされてきた。
一方で、全く手つかずとなっている問題がある。本来は、この地震シミュレーションでは、上記の i)ii)iii)
を一括して解く必要があり、それは 105～6×105～6×104～5 m の領域で、構造物の分解能である最小 10−1m までを考慮した非線形動的問題を解くこととなる。複雑かつ細かい構造物の幾何形状をモデル化し、その表面に沿って応力フリーの境界条件を満足し、応力まで評価し得る手法として、非構造四面体二次要素を用いた陰的時間積分による 3 次元動的非線形有限要素解析での解析を考える。なお、超複雑な幾何形
状を忠実にソリッド要素によりモデル化すると、局所的に小さな要素が生成され、陽的な時間積分では安定性を確保するために非常に小さな時間刻みを使わなければならず、計算量が膨大となるため、安定性に優れる陰解法で解くことも必要となる。結果的にこの問題は、収束性が悪い自由度 1011～13・103～4 時- 50 -

間ステップの以下の陰的時間積分（Newmark-β 法(β=1/4,δ=1/2)）による 3 次元動的非線形有限要素解析を
倍精度で反復型 solver により解くこととなる。
(

42
4
𝑀 + 𝐶 𝑛 + 𝐾 𝑛 ) 𝛿𝑢𝑛 = 𝐹 𝑛 − 𝑄 𝑛−1 + 𝐶 𝑛 𝑣 𝑛−1 + 𝑀 (𝑎𝑛−1 + 𝑣 𝑛−1 ). (2.1.1)
2
𝑑𝑡
𝑑𝑡
𝑑𝑡

ここで、𝛿𝑢、𝑢、𝑣、𝑎、𝐹 は、各々、インクリメント変位、変位、速度、加速度、外力ベクトルである。
また、𝑀、𝐶、𝐾は、整合質量、減衰、剛性マトリクスである。𝑑𝑡、𝑛は時間刻み、タイムステップ数を表している。なお、𝐶 として Rayleigh 減衰マトリクスを用い、その要素減衰マトリクス𝐶𝑒𝑛 は整合要素質量マトリクス𝑀𝑒 及び要素剛性マトリクス𝐾𝑒𝑛 を用いて、Rayleigh 減衰として与える。式(2.1.1) を解いて得られた𝛿𝑢𝑛 により時間ステップ毎に剛性マトリクスを非線形性を考慮して構築し、 𝑄 𝑛 ⇐ 𝑄 𝑛−1 +𝐾 𝑛 𝛿𝑢𝑛 、𝑢𝑛 ⇐ 𝑢𝑛−1 + 𝛿𝑢𝑛 、𝑣 𝑛 ⇐ −𝑣 𝑛−1 + 2 𝑑𝑡 𝛿𝑢𝑛 、𝑎𝑛 ⇐ −𝑎𝑛−1 − 4 𝑑𝑡 𝑣 𝑛−1 + 4 𝑑𝑡 2 𝛿𝑢𝑛 として更新しつ
つ非線形時間発展を解くこととなる。収束性が悪い自由度 1011～13・103～4 時間ステップの式(2.1.1) の求解は、問題の難しさとして上記 i)ii)iii)の state of art を遥かに上回る ultimate earthquake simulation と位置づけられる。このシミュレーションは以前から様々に試みられてきたが、計算量が膨大過ぎるため、i)ii)iii)
を分割して何とか解こうという試みが多数なされている。例えば、特異摂動に基づく階層型解析による[9]とこの流れを汲んだ[6]は、京コンピュータ全系・全力により解析し SC15 の Gordon Bell Prize の finalistにもなっているが、それでも、i)ii)iii)を分割して解くところにようやくたどり着いているのが現状である。
当然のことながら、i)ii)iii)はリンクしており、原理的に切り分けることは難しく、非線形動的問題なので、このリンクを切る行為がどこでどのような影響を与えるかは分からないため、一括解析が望まれる。一方で、上述のように i)ii)iii)を一括して解く ultimate earthquake simulation は計算量が膨大すぎる未達の問題であり、state of art 自体がない状況にある。本研究では、この ultimate earthquake simulation の実現を目指
し、地震シミュレーション分野のフロンティアを切り拓く。
様々なステークホルダーを含む都市・社会の震災軽減では解析結果の信頼性担保が必須となることから、モノづくり等の多くの分野で V & V [13]等による品質保証実績が多数ある非構造四面体二次要素を用いた陰的時間積分による 3 次元動的非線形有限要素解析を対象として solver を開発している。逆に言えば、本 solver 開発の波及効果はこのような多くの分野にも及ぶと期待される。このような幅広な分野で使われるためには、幅広いアーキテクチャにて性能が得られ・他の研究に一般化可能なアルゴリズムであることが重要である。そのため、解の局所領域でのスパースさや解空間のマルチスケールさを活用しつつ高次モードを効率的に推定する 152K-computer-nodes（609,408 MPI processes×12 OpenMP threads =7,312,896 parallel computation)まで良好な性能を示す solver を、なるべく一般的なアルゴリズムを活用しつつ構築することを心がけるとともに、他のアプリでもボトルネックとなる data recurrence を含む計算リッ
チなカーネルの一種である一部の特殊なコアカーネルについては、他のアプリでも参考となるような解決策を、多くのアーキテクチャで高効率な実装が可能なように解決策を示すこととした。また、大規模解析実行中に生成される大規模データから知見を抽出（data learning）し、解析効率向上を目指す様々な方法が提案されているが、良好な超並列効率で良好な結果を得ることは案外難しい。そのため、本研究では、良好な超並列効率を実現しつつ、time to solution を改善し、かつ、他の時間発展問題への導入も容易な data learning 手法も開発し、本 solver に組み込んでいる。以上の詳細を説明するとともに、多様なアーキテクチャにおける豊富な性能計測により、本研究の有効性を以下で示す。

- 51 -

①-(b) 既往研究との比較本研究でターゲットとしている前例のない、ultimate earthquake simulation である、256,000 × 205,000 ×
100,000 m の領域サイズ・構造物の分解能である最小 0.0625m という断層構造物系の非線形動的問題の解
析に、提案 solver である AEGIS (Advanced Earthquake General Integrated Simulator)@富岳で成功した（図2.1.2）。

図 2.1.2 直下型地震に対する超高詳細断層構造物系（関東平野の地盤・地殻と高度に集積されたターミナ
ル駅周辺の都市からなるモデル）を用いた低次非構造 3 次元非線形有限要素解析。モデル全体の領域サ
イズは 256,000 × 205,000 × 100,000 m であり、離散化幅が約 0.0625~8 m で変化する、324,006,449,076 自由度・四面体二次要素数 80,540,873,752 という超大規模問題である。図右半分は都市部の拡大図であり、この領域には地殻の上の柔らかい２層の地盤と 1,060 の建物と 645 のテナントからなる地下街が含まれている。図中央上段は構造物の拡大図であり、最小離散化幅が約 0.0625 m という非常に細かいメッシュ
で構造物の複雑な幾何形状が詳細にモデル化されている。
このような問題を解けると期待された solver として SC14 の Gordon Bell Finalist の低精度演算を活用する京コンピュータ用に開発された GAMERAK [6]があり、これを「富岳」用にチューニングした
GAMERAFugaku[10]、及び、これを高速化した IRISFugaku[11]で解析することを当初計画していた。しかし、今回の問題が巨大すぎかつメッシュサイズの変化が大きすぎるため、as is のままでは反復過程における低精度演算部分の影響により不安定化してしまうことが分かった。様々検討したところ、今回のような
従来にはない巨大な問題でも安定的に良好な超並列効率で解析できる SoA solver は GAMERAK の簡易版としても位置づけられる「京」上で開発された超並列計算に適した PCGEK だった。PCGEK の超並列計算
効率は非常に良いため、これに対して反復過程における低精度演算部分の影響が小さくなるような前処- 52 -

理及び data learning を導入し、「富岳」に最適化することで超並列計算に適した高速な solver AEGIS を開発した。
AEGIS は、「富岳」上での性能計測問題において、SoA である PCGEK と比較して time to solution で 25.45
倍の格段の高速化が達成されている。また、性能計測問題で 9.391 billion DOF@578 computer nodes（実効性能 8.14%）から 1.201 trillion DOF@73,984 computer nodes（実効性能 7.51%）の weak scaling で 93.7%という高い並列性能を示している。また、上記の性能が出にくい現実的な 3,240 億 DOF の適用例問題で、最速の time to solution を目指して 152,352 computer nodes（609,408MPI processes × 12 OpenMP threads =7,312,896 parallel computation) を用いた解析でも、計算ノードあたりの DOF が少なめであるにも関わらず、実効性能 4.66% (21.8 PFLOPS) という良好な性能を実現している。一方で、「富岳」の全系運用は年間で数日程度と予定されているため、多様な地震シナリオを解析可能とするには、「富岳」の通常運用（最大55,296 ノード・24 時間）で適用例が解けることがより重要である。そのため、「富岳」の 52,992 計算ノード（211,968 MPI processes × 12 OpenMP threads = 2,543,616 parallel computation）を使用して適用例の解析（9,550 時間ステップ，dt=0.005s）も行ったところ、4 時間 36 分という短時間で解析することが出来た。
なお、SoA である PCGEK を「富岳」で用いてこの適用例を計算しようとすると 95.8 時間かかる見込み
であり、解析実行が難しいことがわかる。実アプリを使った検討実行の観点からも、AEGIS の有効性が際立つ結果となっている。AEGIS、SoA である PCGE、及び上記の関連する solver について、適用例及び適用例を模した性能計測問題による詳細な性能評価は、Section①-(e)にて示す。
①-(c) 開発手法
幅広い分野で使われている 3 次元有限要素解析の式(2.1.1) を超並列で高速に解くための Algorithm 1を開発した。具体的には、アーキテクチャ依存性の低い一般的な演算をベースとした一般的な implicit
solver である共役勾配法を拡張し、解の局所領域でのスパースさを活用する高性能な initial solution
predictor 及び解空間のマルチスケールさ及び低精度演算を活用する preconditioner を組み込むことで、収束性を悪化させる原因となる高次モードの誤差を効率的に予測・減少させ、低精度演算の影響を受けにくくしつつ、massive parallel computing を可能とすることで上記を実現している。また、このアルゴリズ
ム を 実 現 す る 上 で ア ー キ テ ク チ ャ 依 存 性 が 生 じ る 一 部 の 特 殊 な カ ー ネ ル （ K、 C 、 M に 関 す る
Element-by-Element (EBE) kernel）が必要となり、低精度演算下での高速化が肝となるが、この各アーキテ
クチャで高効率な実装が可能なことを示すことで、Algorithm 1 がロバストに多くのシステムで高効率が期待できることを示した。なお、EBE は、一般のアプリに頻出する data recurrence を含む計算リッチでボトルネックとなるカーネルの一種であるため、他のアプリにも有用となるように解決の一方策を示している。以下に各詳細を示す。

- 53 -

①-(c) A. Initial solution estimator
Implicit solver の反復回数の低減に直結する、より良い initial solution を推定する predictor へのニーズは高い。ポイントは超並列下での解析コストとその推定解の性能のバランスである。ここでは、dynamic
mode decomposition (DMD)[12] の考えを拡張した、𝑛 − 1 step までの計算済みの結果を学習し、𝑛 step 目の初期解を推定する、massive parallel computation に適し、かつ、低解析コストの initial solution estimatorを構築する(Algorithm 1 の line 3 で用いる)。
DMD は、解の時系列データを用いて近似逆行列を介してその時間発展を表すオペレータを推定する手法であり、このオペレータを用いることで、現在の解から次の解を予測できると期待される。一方で、
transit な挙動を含む多数のモードを DMD により推定するには解析コストが多くかかるとともに、trendを含む場合には stacking が必要になる場合もあり、さらに解析コストがかかる。
以上を踏まえ、trend を含む低次モード解は解析コストの小さい他手法で推定し、残りの高次モードを局所的な時空間において推定することを超並列計算に適した形で実現する。具体的には、解析領域を小𝑛領域に分割し、各領域において Adams-Bashforth method により、trend を含む低次モード解𝛿𝑢𝑎𝑑𝑎𝑚を直近
𝑑𝑡

4 ステップの速度を用いて24 (−9𝑣 𝑛−4 + 37𝑣 𝑛−3 − 59𝑣 𝑛−2 + 55𝑣 𝑛−1 )により推定する。真の解𝛿𝑢𝑛 との差𝑥 𝑛- 54 -

𝑛を𝛿𝑢𝑛 − 𝛿𝑢𝑎𝑑𝑎𝑚と定義し、過去の𝑠 + 1 step 分の解析結果を用いて、𝑋 𝑛−1 = [𝑥 𝑛−1 𝑥 𝑛−2 … 𝑥 𝑛−𝑠 ]を定義し、この行列を用いて、時間発展を表すオペレータ𝐴を𝑋 𝑛−1 = 𝐴𝑥 𝑛−2として定義する。Trend を含む低次モー𝑛𝑛−1𝑛) + 𝛿𝑢𝑎𝑑𝑎𝑚
ドを含み、かつ、高次モードも含む初期解𝛿𝑢𝑖𝑛𝑖
は、𝐴(𝛿𝑢𝑛−1 − 𝛿𝑢𝑎𝑑𝑎𝑚
として計算される。な𝑛−1)の計算は Modified Gram-Schmidt orthonormalization により行う。
お、DMD とは異なり𝐴(𝛿𝑢𝑛−1 − 𝛿𝑢𝑎𝑑𝑎𝑚

以上のように、小領域毎に推定するので、その挙動を予測するために必要なモード数（すなわち、学習に必要なデータ数）を軽減することが出来、かつ、大域的な計算を避けることで高効率な超並列計算が可能としている。
𝑛
学習データが十分にある場合、例えば、性能計測問題では𝛿𝑢𝑖𝑛𝑖
は相対誤差で 10−5 程度の精度で初期解𝑛を推定出来ている（𝛿𝑢𝑎𝑑𝑎𝑚
は 10−3 程度の精度）。一方で、学習データが十分でない場合には、当然のこと
𝑛𝑛ながら、推定精度は悪化するため、そのような場合には𝛿𝑢𝑖𝑛𝑖として、𝛿𝑢𝑎𝑑𝑎𝑚を用いるように自動的に切り替わるようにしている。
①-(c) B. Preconditioner for r
ベースとなる共役勾配法として、
Algorithm 1 は、adaptive conjugate gradient の一種である inexact conjugate
gradient [20]としており、adaptive conjugate gradient の特性を活かし、𝐴𝑛 𝑟 = 𝑧を単精度で粗い閾値で CG により解くことを preconditioner とする(Algorithm 1 の line 8 で用いる)。具体的には、𝐴𝑛 𝑟 = 𝑧を解く𝐶𝐺𝑝𝑟𝑒 、
𝐴𝑛 𝑟 = 𝑧 の lower mode を 推 定 す る こ と を 目 的 と し 𝐴𝑛 𝑟 = 𝑧 と 等 価 な ラ フ な 問 題 を 解 く CG を
𝑙𝑚𝑟
ℎ𝑚𝑟
𝐶𝐺𝑝𝑟𝑒
、𝐴𝑛 𝑟 = 𝑧の一部の高次モードが卓越する領域のみを重点的に解き higher mode refine を行う𝐶𝐺𝑝𝑟𝑒
𝑙𝑚𝑟
𝑙𝑚𝑟
ℎ𝑚𝑟
を用いて、𝐶𝐺𝑝𝑟𝑒 、𝐶𝐺𝑝𝑟𝑒、𝐶𝐺𝑝𝑟𝑒 、𝐶𝐺𝑝𝑟𝑒、𝐶𝐺𝑝𝑟𝑒の順で解く修正 W サイクルマルチグリッドにより解く。

これは、Preconditioner であるため、粗く解いて良く、低次モードは高次モードに比べラフな取り扱いをしても大きな問題は生じない。一方で、高次モードが精度よく考慮されるように解かないと収束性が悪化するので、今回のような収束性の悪い高次モードが含まれるような場合には解析コストバランスを考
𝑙𝑚𝑟
慮した適切な設計が必要となるためである。本研究では、𝐶𝐺𝑝𝑟𝑒
は𝐴𝑛 𝑟 = 𝑧を四面体一次要素としてブロックヤコビ前処理を用いて単精度で解く。巨大な DOF ではあるものの、低次モードを扱っていることとから、そこまで精度の心配はしなくて良く、かつ、𝐶𝐺𝑝𝑟𝑒 よりは自由度も少ないため、対象行列は CRS
𝑙𝑚𝑟
𝑙𝑚𝑟
で保存して計算する（𝐶𝐺𝑝𝑟𝑒
は CRS での一般的な CG である）。𝐶𝐺𝑝𝑟𝑒の反復回数は多いため、CRS 構築

コストを考慮しても高速に計算可能となる。一方で、高速に計算しつつも高次モードの精度劣化を抑えℎ𝑚𝑟ることが必須の𝐶𝐺𝑝𝑟𝑒 及び𝐶𝐺𝑝𝑟𝑒
は𝐴𝑛 に関する計算に注意が必要となる。Initial solution estimator により高次モードがすでに精度よく推定されているため、ブロックヤコビ前処理を用いて解くこととしているとℎ𝑚𝑟
はいえ𝐶𝐺𝑝𝑟𝑒 及び 𝐶𝐺𝑝𝑟𝑒の反復回数は十分少なくなっていること、また、自由度が大きいため CRS の構

築コストが大きくなることから、局所的な演算を積み上げるため精度劣化抑制も期待される高速な EBEℎ𝑚𝑟
による𝐴𝑛 が必要となる（なお、𝐶𝐺𝑝𝑟𝑒 及び𝐶𝐺𝑝𝑟𝑒
においても𝐴𝑛 に関する以外は一般的な CG である）。しかし、EBE は data recurrence を含む計算リッチなカーネルかつアーキテクチャ依存性の高いカーネルであるため特別な配慮が必要となってしまうが、以下に示すように汎用的な実装法があるので問題ない。
①-(c) C. Data recurrence を含む計算リッチなカーネルの実装
全体行列を構築し行列ベクトル積計算時に全体行列をメモリから読みだす通常の行列ベクトル積計算
に対し、メモリフリーの行列ベクトル積計算手法である EBE 法では、行列ベクトル積計算の都度要素行- 55 -

列𝐴𝑛𝑒 を構築し要素右辺ベクトルに掛け合わせた上で(𝑓𝑒𝑛 ⇐ 𝐴𝑛𝑒 𝑁𝑒𝑇 𝛿𝑢𝑛 ) 左辺ベクトルに結果を足しこむ過程をとる(𝑓 𝑛 ⇐ ∑𝑒 𝑁𝑒 𝑓𝑒𝑛 )。ここで、𝐴𝑛𝑒 は節点位置と物性情報から計算される要素マトリクスであり、𝑁𝑒は要素ローカル節点番号をグローバル節点番号にマッピングする行列である。EBE 法においては要素毎にローカルに変数をスケーリングできるため、低精度演算利用時において変数のオーバーフロー・アンダ
ーフローをコントロールしやすく、また、要素計算時に近い値の差分をとり得る変数においてはあらか
じめ shifting を実施することで精度劣化を抑えることが可能である（例えば、要素剛性マトリクスにおいては右辺ベクトルの平行移動成分とその他のモードに分けることで桁落ちを抑えることが可能である）。
さらに、EBE 法においては特に高い精度が必要な一部の演算(e.g., 座標値計算) のみに高精度演算を選択的に適用できるため、高精度演算を実施するためには全ての演算を倍精度化する必要のある全体行列格
納形式の行列ベクトル積と比べて柔軟に精度の調整が可能である。その一方で、EBE 法は全体行列格納
形式の行列ベクトル積と比べて全演算量が多く高効率計算が求められるが、data recurrence を含むランダ
ム足しこみが含まれるため straight forward に並列計算することができない。そこで、本研究ではメニー
コア・wide SIMD CPU と高速な atomic 演算が可能な GPU を対象に各計算機アーキテクチャの特性に合わせてカーネルのローカルなアルゴリズムを設計することで汎用的な実装を実現する。
以下、図 2.1.3 の EBE 法の serial code を用いて詳細を説明する。line 18–20 における f のインデックス(cny1、cny2、...、cny4)が不規則性・再帰性を持つため、ベクトル f への値の足しこみ時に data recurrenceが生じ、ie ループでの並列化が阻害される。図では省略しているが line 14–16 の計算は実際には数百行からなるループ内に閉じた計算となっており、ここが並列計算できない場合はカーネルの性能が著しく劣化する。そこでまず、CPU における ie ループのマルチスレッド並列化を考える。一番単純な方法は、足しこみ先 f をスレッド数だけ確保し(f1、f2、...、fnthread)、全要素計算後にそれらのスレッド毎の結果を全体ベクトル f に足しこむ方法である(f =Σnthread fi)。実装は容易であるが、コア数が多いメニーコア環境においてはメモリ利用量が多く最後のスレッドローカルな結果の集約に多大なコストが生じることとなる。
そこで、[10]にならってキャッシュ効率を保ちつつマルチコア間でデータ競合が生じないようなスレッド
分割方法を開発した(図 2.1.4b)。この方法は図 2.1.4a に示す通常のカラーリングのように各カラー内の全ての要素間で節点の重なりが無いように領域分割するのではなく、各スレッド内が受け持つ領域内ではrecurrence を許すため、キャッシュの再利用率を高めることが可能である。併せて、要素ループボディをデータ読み込み＋演算部と、データ足しこみの 2 つのループに分割することでデータ読み込み＋演算部
において SIMD 演算を使えるようにした(図 2.1.5)。ループブロックの大きさ NL を SIMD 長と同じ値に設定することで、テンポラリバッファ fe11、fe21、...fe34 をレジスタ・L1 キャッシュに載せたまま計算
できることとなる(i.e., 512-bit SIMD CPU で単精度演算の場合は NL=16、256-bit SIMD CPU で単精度演算の場合は NL=8)。このようにすることでメニーコア・wide SIMD 演算器を用いることが可能となっており、幅広い CPU アーキテクチャで性能向上が期待できる。アーキテクチャのさらに詳細な特性にあわせて高度化することも可能であり、例えば、「富岳」の A64FX CPU のレジスタ数・out-of-order 資源・演算/デー
タアクセスレイテンシにあわせてループのさらなる分割や、分割されたループ間での演算スケジューリ
ングの制御により高速化を達成している [18]。
GPU においては超多数のスレッドを同時に立ち上げることでデータアクセス・演算レイテンシを隠蔽することが性能達成のために必須となるが、上記の CPU 版と同じスレッド分割・ループ分割アプローチ
ではスレッド数が多すぎて性能が十分に発揮されないこととなる。その一方で、近年の NVIDIA GPU に- 56 -

は高速な hardware accelerated atomic add 命令が実装されているため、これを活用したカーネルを実装した(図 2.1.6)。ここでは、data recurrence を含む f への足しこみ部のみ atomic add を用いるようにしている(line20–25)。この atomic add は分散配置された L2 キャッシュに隣接した演算器で実行可能なため、コアまでのデータ移動の必要なく高速に演算可能である。なお、全節点の情報を直接 global memory に atomic add
せずに、複数のスレッドで共有される高速メモリ(e.g., shared memory)上で f の値を集約し、その後集約した値を global memory に atomic add することでさらなる高速化も可能である。Atomic add がサポートされていない GPU においては上記のローカルな値の縮約と図 2.1.4b の領域分割法を組み合わせることで高速に EBE 計算が可能であり、アーキテクチャに応じ上記のいずれかを選択することで効率よく EBE カーネ
ルが演算できると期待される。

図 2.1.3 Serial EBE kernel structure for linear tetrahedral element. The code structure is the same for thetargeted second order tetrahedral element but with 10 × 3 = 30 components of ue and fe per element.

- 57 -

図 2.1.4 Coloring and thread partitioning methods for EBE computation- 58 -

図 2.1.5 Tuned EBE kernel for many-core wide SIMD architecture- 59 -

図 2.1.6 Tuned EBE kernel for GPU with hardware accelerated atomics①-(c) D. Massive parallel computation
超大規模問題の求解時における性能低下要因となる、計算量の増大・ロードインバランスの発生・通信コストの増大、をそれぞれ抑える工夫を実施した。
計算量増大の抑制：本研究での求解アルゴリズムはポイントワイズのブロックヤコビ前処理ベースのCG 法で構築されているため、各反復における計算量は問題規模(要素数) に比例する。従って、問題規模によらず CG 法の反復数を一定に保つことができれば、大規模問題まで weak scaling が達成できると期待される。低精度演算による解の精度劣化が生じると高次モードをとらえきれなくなり、反復数の増加が懸念されるため、EBE 法における行列ベクトル積の精度劣化抑制策に加え、MPI_Allreduce 時の桁落ち回避のための入力ベクトルのスケーリングを実施した。また、前処理の w-cycle における各前処理
CG の反復数を抑えることで倍精度演算の同期回数を増やし解が発散するのを防いでいる。
𝑙𝑚𝑟ℎ𝑚𝑟
ロードインバランスの抑制：𝐶𝐺𝑝𝑟𝑒 ・𝐶𝐺𝑝𝑟𝑒
は対象領域全体をカバーするのに対し、𝐶𝐺𝑝𝑟𝑒
は𝐶𝐺𝑝𝑟𝑒 のう
𝑙𝑚𝑟
ち収束性の悪い一部の領域を選択的に解くこととなるため、MPI 計算時においては𝐶𝐺𝑝𝑟𝑒・𝐶𝐺𝑝𝑟𝑒
、及び、
ℎ𝑚𝑟
𝐶𝐺𝑝𝑟𝑒
の双方でロードバランスが取れるように領域分割をする必要がある。そこで、メッシュの接続情報ℎ𝑚𝑟ℎ𝑚𝑟と要素物性情報からあらかじめモデルのうち収束性の悪い部分を特定し𝐶𝐺𝑝𝑟𝑒
部を選択した上で、𝐶𝐺𝑝𝑟𝑒- 60 -

部とそれ以外の 2 領域においてそれぞれグラフ分割ソフトウェア[19] を用いて領域分割を実施すること
で、超並列環境においても通信量を抑えつつロードバランスを保っている。
通信コストの抑制：超並列環境においては、大域通信(MPI_Allreduce) に加え、通信ホップ数が増える
ために隣接通信のレイテンシも増加する傾向にあり、また、他ジョブのノイズの影響を受けやすくなる傾向にある。過去の求解データの学習による初期解の高精度推定、及び、coarse-grid による効率的な長波長成分の求解により、通信の総呼び出し回数を削減し、レイテンシ・ノイズの影響を抑える工夫をしている。また、「富岳」においては複数ノードの IO を処理する機能を一部の計算ノードが持ち合わせているため、IO ノードを専有利用するノードマッピングでジョブ実行することで他ジョブのノイズ混入を防いでいる。
①-(d) 性能計測問題・性能計測環境
Section ①-(b)で示した適用例を模した、断層から都市の構造物までを含む有限要素法モデルを生成し、非線形地震波動計算にかかる時間を計測した。ここでは、東京周辺の実際の地殻・地盤・構造物のデジ
タルデータを元に作成した 65,500×65,500×6,000m の地殻・都市統合モデルを用いる（図 2.1.7）。都市部は構造物の複雑形状を解像し、軟弱地盤が非線形化しても数値解の収束が担保されるよう最小要素サイズ
0.25m、地殻部は地殻構造を反映した数値解の収束が担保できる最小要素サイズ 32m でモデル化している。
1995 年の兵庫県南部地震において JMA 神戸で観測された地震動[21]の振幅を 1/10 にしたものを時間ス
テップ幅は dt=0.005s でモデル底面に入力する。図 2.1.7 に示すように対象のモデルは幅広い物性と要素
サイズを持つため、従来の地殻のみ・都市のみの地震シミュレーションと比べて性質が悪く、反復法で効率的に求解するのはチャレンジングな問題となっている。Weak scaling を測るため、上記モデルを x, y方向に鏡面を取って接続することで、類似の収束特性を持った大規模モデルを作成する(表 2.1.1)。地震動
が都市部に伝わり反復解法の収束性が悪化しはじめる 401 から 600 時間ステップにおける反復数・実行時間を計測する。

- 61 -

図 2.1.7 Performance measurement setting, created based on digital data of the actual crust, ground, andstructures.
表 2.1.1 Configurations of the models used to measure performance. Model set A is a weak scaling model set
while model set B is a strong scaling model set.

性能計測問題においては提案手法及び SoA である PCGE の性能を比較する。PCGE は Algorithm 1 における line 3 の初期解推定に Adams-Bashforth method を用い、line 8 の preconditioner に倍精度の 3×3 ブロッ
クヤコビ前処理を用いた solver に相当する。主要演算となる行列ベクトル積が EBE 法にて実行されるためロードバランスに優れ、超並列環境において幅広い大規模問題に対してロバストに高速求解できる手法となっている。提案手法における predictor の領域分割においては、一領域当たりの四面体二次節点数
が約 20,000 節点となるようグラフ分割ソフトウェアを用いて各 MPI 領域をさらに分割した。なお、非線形構成則の評価にかかる時間は対象問題によって異なるため、他問題においても汎用的に比べることが可能な初期解推定を含む solver 部に要した時間・演算性能を比較する。開発手法の超並列時の実行性能- 62 -

を「富岳」上で計測する。加えて、手法が幅広いアーキテクチャにおいて有効であることを示すため、Xeon (Cascade Lake) CPU ベースの Oakbridge-CX システムでも性能を比較する。また、カーネルレベル
では、NVIDIA A100 GPU (PCI-E 版)、V100 GPU (NVLink 版)、AMD EPYC (Rome) CPU においても性能を計測する。
「富岳」は 48 計算コアを有する Arm v8.2-A SVE ベースの A64FX CPU [15] 1 個を搭載する計算ノード
158,976 台からなるシステムである[14]。各コアの SIMD 演算器は 512 ビットであり、各ノードのピーク
性能は倍精度で 32 x 48 コア x 2GHz = 3,072 FP64 GFLOPS となる(2.2 GHz のクロックブーストで 3,379
FP64 GFLOPS)。FP32 は FP64 のハードウェアピークの倍となる。各計算ノードは HBM2 メモリ 32 GB を搭載し（ピークメモリバンド幅は 1024 GB/s）、各計算ノードは 6 次元トーラスネットワーク(TofuD [16])でつながれている。A64FX CPU は 12 コア(core memory group; CMG)毎にまとまった NUMA アーキテクチ
ャとなるため、本研究では、4 つの CMG においてそれぞれ MPI プロセスを 1 個立ち上げ、そこからそれ
ぞれ 12 OpenMP スレッドを立ち上げる方法で並列計算をした。
Oakbridge-CX (OBCX) [17]は東京大学情報基盤センターに設置されている Xeon CPU server ベースのスーパーコンピュータシステムである。1,368 台ある計算ノードはそれぞれ 28 コア Xeon Platinum 8280
（Cascade Lake）CPU を 2 台、192 GB DDR4 memory (281.6 GB/s bandwidth)を搭載する。各 CPU は 512-bit
SIMD units を搭載しているため、ノード当たりピーク性能は 4.84 TFLOPS となる。各計算ノードは
full-bisection fat tree topology をもつ Intel Omni-Path Architecture でつながれている。
①-(e) 性能計測結果
本手法が多種の計算機アーキテクチャにおいて性能が発揮されることを示すため、現時点で利用できる主要な高性能計算アーキテクチャにおける性能を計測する。その後、本手法の超大規模問題における性能を確認するため、超多数計算ノードを利用可能な「富岳」において weak scaling, 及び、大規模問題における strong scaling を計測する。最後に、適用例での性能計測を通して実問題での性能を確認する。
①-(e) A. 主要な高性能計算アーキテクチャにおける性能
まずはアーキテクチャ依存性の高い EBE カーネルにおける性能を代表的な CPU アーキテクチャ、及び、GPU アーキテクチャにおいて計測する。表 2.1.2 に「富岳」の A64FX CPU、Xeon CPU、EPYC CPU、
A100 GPU、及び V100 GPU における EBE カーネルの性能を示す。ここでは、model A-1 のプロセス番号
0-3 のデータを用いて計測を実施している。SIMD 幅が 2 と狭くコア数が 8 であった「京」用に SC14 のGBF で開発された as is の EBE カーネルに比べ、キャッシュ効率と wide SIMD を考慮して開発した EBE
カーネルでは実行性能が顕著に伸び、「富岳」の A64FX では倍精度ピーク比で 12.3%、Xeon と EPYC CPU
では 21-22%の倍精度ピーク比となった。また、CPU とは異なるアプローチで EBE を計算する GPU 版
カーネルを A100 及び V100 GPU で実行した際の倍精度ピーク比はそれぞれ 22%、19%となった。このように、計算アーキテクチャに合わせたカーネルアルゴリズムを採用することで幅広いアーキテクチャにおいて高い性能が達成されていることがわかる。なお、一見 A64FX の性能が悪いように見えるが、これは電力性能の向上のために out-of-order リソースなどを絞り、また、メモリバンド幅を十分にとった設計となっているためであり、後述のようにメモリバンド幅ネックのカーネルが含まれるフルアプリではピーク性能比で優れた性能となっている。

- 63 -

表 2.1.2 多種のアーキテクチャにおける EBE カーネル(倍精度四面体二次要素) の性能。それぞれ 4 プ
ロセスで計測した場合の実行時間(ms) とハードウェアピークに対する FLOPS の比率を示している。

次に、開発手法、及び、SoA である PCGEK の solver 全体の性能を「富岳」及び Xeon ベースの OBCXにおいて計測する(図 2.1.8)。ここでは、「富岳」578 ノード（倍精度ピーク 1.77 PFLOPS、メモリバンド
幅 591 TB/s）、及び、OBCX 256 ノード(1.24 PFLOPS, 72.1 TB/s) を用いて model A-1 を解析した。後述のように演算数が PCGEK 比で 1/8.16 に削減され、さらに、アーキテクチャに適した EBE カーネルを用いることで高い演算効率（「富岳」において倍精度ピーク比で 8.14%、OBCX においては 7.12%）が得られているため、PCGEK 比で「富岳」で 26.0 倍、OBCX で 18.9 倍という高い高速化率が得られている。以後、性能向上要因の詳細を分析する。401～600 時間ステップを求解する際の PCGEK solver の総反復数はℎ𝑚𝑟
𝑙𝑚𝑟
62,345 回だったのに対し、AEGIS の合計反復数は𝐶𝐺𝑝𝑟𝑒 4,080 回、𝐶𝐺𝑝𝑟𝑒
8,160 回、𝐶𝐺𝑝𝑟𝑒
30,760 回の合計 43,000 回であった。PCGEK で用いられている Adams-Bashforth method による初期解の推定精度（相対誤差 10−3 程度) が、提案手法の predictor により大幅に向上していること（相対誤差 10−5 程度）、また、
𝑙𝑚𝑟
提案手法においては lower frequency mode を multi-grid における coarse grid (𝐶𝐺𝑝𝑟𝑒）で効率よく求解できていることが、反復数削減につながっていると考えられる。PCGEK の CG 反復、及び、提案手法におけℎ𝑚𝑟
𝑙𝑚𝑟
る𝐶𝐺𝑝𝑟𝑒 、𝐶𝐺𝑝𝑟𝑒、𝐶𝐺𝑝𝑟𝑒の各反復の自由度の比が 1 : 1 : 0.087 : 0.125 となることを考慮すると、反復数の比率以上に演算数が減少することとなり、実際、演算数は 1/8.16 となった。図 2.1.9 に示す実行時間内訳
𝑙𝑚𝑟ℎ𝑚𝑟
を見ると狙い通り𝐶𝐺𝑝𝑟𝑒と𝐶𝐺𝑝𝑟𝑒における反復あたりの実行時間が𝐶𝐺𝑝𝑟𝑒 比で抑えられていること、また、

predictor にかかる時間は全体の 1/10 程度に抑えられていることがわかり、これらの反復数削減・反復当たりの演算数削減・演算効率の向上、の相乗効果が「富岳」
・OBCX 双方における高速化につながってい
ることがわかる。

- 64 -

図 2.1.8 Model A-1 における「富岳」
・ Oakbridge-CX（OBCX）での solver+predictor の実行時間比較。
カッコ内に倍精度ピーク性能比を示す。

図 2.1.9 Breakdown of elapsed time of AEGIS on Model A-1. Note that the predictor is data access bandwidthbound, and thus relatively high performance is achieved on Fugaku with higher memory bandwidth.
①-(e) B. 大規模問題における性能
次に本手法の大規模問題における性能を、現時点で最大規模の計算環境である「富岳」において計測する（図 2.1.10）。まず、提案手法は、他手法と比較して、良好な time-to-solution 及びスケーラビリティを示していることが分かる。SoA とされる PCGEK と比較して、A-3 においても 9,318/366 ＝ 25.45 倍の高速化を達成している。PCGEK は時間がかかりすぎ大規模問題においては計算資源の都合上計測できな
いため、本研究において開発した EBE カーネルを組み込んで SoA を高速化した PCGEFugaku を開発し、より大規模問題を解いたが、PCGEK から PCGEFugaku で A-3 において 9,318/1,718 ＝5.42 倍の高速化に留まっており、提案手法に比べると相当遅く、提案手法の高速さが際立つ結果となっている。一方で、当初
SoA として期待していた GAMERAK を本研究において開発した EBE カーネルを組み込んで開発した
GAMERAFugaku にて大規模問題を解いてみたが、問題規模・計算ノード数の拡大に伴い反復数が増大し、- 65 -

36,992 ノード以上では solver が収束しなくなるなど、本研究でターゲットとする超大規模問題の求解自体が不能となっている。これは、大規模問題求解時の低精度演算のケアの問題であり、実際、GAMERA の前処理を全て倍精度で計算することで大規模問題における反復数の増加は抑えられることがわかったた
め、GAMERA の手法自体に問題があったわけでなく、今回対象とするかなり複雑で特性が大きく変化する問題においては低精度演算による前処理の精度が不足していたことがわかった。なお、GAMERAK 及び GAMERAFugaku は小規模問題であっても遅く、大規模問題は解けないし遅いということで、提案手法に比べると劣ってしまっている。その一方で、低精度演算による誤差拡大を抑える本手法においては全工程を倍精度で計算する PCGEFugaku と同様に高並列まで反復数がフラットに維持され、大規模モデルにおいても高い weak scaling 性能が達成されていることがわかる。提案手法は PCGE 比で総反復数とそれに伴う MPI_Allreduce や隣接通信の回数が減っているため、提案手法の PCGEFugaku 比の実行時間は最小のmodel A-1 での 4.73 倍速から比較できる最大の model A-7 で 5.04 倍速に拡大している。Model A-1 とmodel A-8 での weak scaling は 93.7%であり、model A-8 でのピーク性能は 17.07 PFLOPS (ハードウェアピ
ークの 7.51%)であった。また、図 2.1.11 に示す性能計測問題における多数ノードへの計測結果が示す通り本手法は高い strong scaling 性能が得られており、ノード数に対して問題規模が小さい場合でも PCGE
比で高い高速化率が実現されている。

図 2.1.10 Size up scalability for model set A.

- 66 -

図 2.1.11 Speedup scalability for model set B. Note that only problem sizes up to B-3 is used for actual use
purposes and thus B-4 is only for a reference aimed for shortest time-to-solution.
①-(e) C. 実際のアプリケーションにおける性能
Section ①-(b)にて示した性能が出にくい現実的な 3,240 億 DOF の適用例を例に AEGIS の性能を確認する。まず「富岳」全系で適用例解析の最速の time to solution を狙った場合を考え、152,352 computer nodes（609,408 MPI processes×12 OpenMP threads = 7,312,896 parallel computation) で解析してみたところ、
computer node あたりの DOF が少な目であるにも関わらず実効性能 4.66% (21.8 PFLOPS) という良好な性能を実現していることがわかった。一方で、「富岳」の全系運用は年間で数日程度と予定されているおり、多様な地震シナリオを解析可能とするには、「富岳」の通常運用（最大 55,296 ノード・24 時間）で適用例が解けることがより重要であるため、「富岳」52,992 ノードでの AEGIS の性能も確認する。比較対象として、PCGEK を用いたいが、「富岳」上であっても遅すぎるため、PCGEK に本研究で開発した EBEを導入して高速化した PCGEFugaku を用いる。比較可能な 1,500–5,700 タイムステップでの平均で、提案手法は PCGEFugaku の 4.4 倍速となったことから、SoA である PCGEK との比較がもしできれば、
AEGIS は 23.87
倍速となると推察される（性能計測問題での PCGEFugaku と PCGEK の 9318 s/1718 s = 5.42 倍の差によりス
ケーリングした）。また、AEGIS は 9550 ステップまでの解析を 4 時間 36 分で実行しており、1 タイム
ステップ当たり 1.5 秒と極めて高速に解析出来ている。一方で、SoA である PCGEK を用いると 4.0×23.87= 95.8 時間かかると推定され、SoA では「富岳」であっても今回のような適用例を解析することは難し
いことが分かる。なお、この適用例での AEGIS/PCGEFugaku の高速化率は性能計測問題での高速化率(4.7～5.0 倍)とほぼ同じであり、震源過程が性能計測問題よりもさらに複雑な実用的な適用例問題においても本手法がロバストに高速求解できることが示された。これは、提案手法の predictor による初期解の誤差がほぼ全ての時間ステップにおいて Adams-Bashforth method よりも 1～2 桁低い誤差に抑えられていること、また、低精度演算によるマルチグリッド前処理が効いていることに起因する。

①-(f) まとめと今後の展望- 67 -

152 K-computer-nodes（609,408 MPI processes × 12 OpenMP threads = 7,312,896 parallel computation) まで良好にスケールする implicit solver を開発し、3 次元低次非構造要素を用いた動的非線形有限要素解析による超大規模解析を実現した。
「富岳」上での性能で SoA と比較して 25.45 倍の高速化を達成した提案手法により、SC15 の Gordon Bell Finalist の solver@京ですら階層型解析による近似解析しか出来なかった問
題の直接解析を実現するに至った。なお、SoA@京と提案手法@富岳を比較すると、ハードウェア性能 42倍×手法の性能差 25.4 倍＝1,070 倍の解析能力差となっている。断層から構造物までの一括解析は地震工学の夢の一つであり、それが実現可能になったことを示した本解析は、今後の地震工学のシミュレーシ
ョンにおけるマイルストーンの一つとなると考えられる。また、本手法は収束性の悪い高次モードを含
む時間発展問題に援用可能であり、多くの問題で高速化が期待される。特に、解析途中で生成されるデータを用いた超並列可能な data learning による initial solution estimator は、簡便に他の solver に移植することが可能な汎用性の高いものである。また、他のアプリでも頻出する data recurrence を含む計算リッチなカーネルを高速化する方法を、各アーキテクチャ上で汎化性のあるように示したことは、本 solver の汎用性を示すだけでなく、他のアプリの同様のカーネルの高速化にも寄与すると期待される。

[1] T.H. Jordan, SC15 Invited Talk Spotlight: Societal Impact of Earthquake Simulations at Extreme Scale.
[2] L. Carrington, D. Komatitsch, M. Laurenzano, M. Tikir, D. Michea, N.L. Goff, A. Snavely and J. Tromp,
High-frequency simulations of global seismic wave propagation using SPECFEM3D_GLOBE on 62K
processors, Proceedings of the International Conference on High Performance Computing, Networking, Storage
and Analysis (SC’08). IEEE Press, 2008, pp. 1–11.
[3] A. Heinecke, A. Breuer, S. Rettenberger, M. Bader, A-A. Gabriel, C. Pelties, A. Bode, W. Barth, X-K. Liao, K.
Vaidyanathan, M. Smelyanskiy, and P. Dubey. “Petascale High Order Dynamic Rupture Earthquake Simulations
on Heterogeneous Supercomputers.” Proceedings of the International Conference on High Performance
Computing, Networking, Storage and Analysis (SC’14). IEEE Press, 2014, pp. 3–14.
[4] H. Fu, C. He, B. Chen, Z. Yin, Z. Zhang, W. Zhang, T. Zhang, W. Xue, W. Liu, W. Yin, G. Yang, X. Chen,
15-Pflops Nonlinear Earthquake Simulation on Sunway TaihuLight: Enabling Depiction of Realistic 10 Hz
Scenarios, SC ’17 Proceedings of the International Conference for High Performance Computing, Networking,Storage and Analysis, Article No. 2.
[5] C. Uphoff, S. Rettenberger, M. Bader, E. H. Madden, T. Ulrich, S. Wollherr, and A.-A. Gabriel, “ Extreme scale
multi-physics simulations of the tsunamigenic 2004 sumatra megathrust earthquake, ” SC ’17 Proceedings of theInternational Conference for High Performance Computing, Networking, Storage and Analysis, pp. 1–16, 2017.
[6] T. Ichimura, K. Fujita, S. Tanaka, M. Hori, M. Lalith, Y. Shizawa, and H. Kobayashi. “Physics-based urban
earthquake simulation enhanced by 10.7 BlnDOF x 30 K time-step unstructured FE non-linear seismic wave
simulation.” Proceedings of the International Conference on High Performance Computing, Networking, Storage
and Analysis (SC’14). IEEE Press, 2014, pp. 15–26.
[7] T. Ichimura, K. Fujita, P. E. B. Quinay, L. Maddegedara, M. Hori, S. Tanaka, Y. Shizawa, H. Kobayashi, and K.
Minami. “Implicit Nonlinear Wave Simulation with 1.08T DOF and 0.270T Unstructured Finite Elements toEnhance Comprehensive Earthquake Simulation.” Proceedings of the International Conference on High
Performance Computing, Networking, Storage and Analysis (SC’15). ACM, 2015, p. 4.

- 68 -

[8] T. Ichimura, K. Fujita, T. Yamaguchi, A. Naruse, J. C. Wells, T. C. Schulthess, T. P. Straatsma, C. J. Zimmer, M.
Martinasso, K. Nakajima, M. Hori, L. Maddegedara. “A fast scalable implicit solver for nonlinear time-evolution
earthquake city problem on low-ordered unstructured finite elements with artificial intelligence andtransprecision computing.” Proceedings of the International Conference on High Performance Computing,
Networking, Storage and Analysis (SC’18). ACM, 2018, Article No. 49.
[9] T. Ichimura, M. Hori, P. E. B. Quinay, L. Maddegedara, T. Suzuki and S. Noguchi, Comprehensive numericalanalysis of fault-structure systems – Computation of the large-scale seismic structural response to a given
earthquake scenario –, Earthquake Engineering & Structural Dynamics, 41, pp 795–811, 2012 April, DOI:
10.1002/eqe.1158.
[10] Kohei Fujita, Kentaro Koyama, Kazuo Minami, Hikaru Inoue, Seiya Nishizawa, Miwako Tsuji, Tatsuo Nishiki,Tsuyoshi Ichimura, Muneo Hori and Lalith Maddegedara, High-fidelity nonlinear low-order unstructured
implicit finite-element seismic simulation of important structures by accelerated element-by-element method,
Journal of Computational Science, 2021.
[11] T. Ichimura, K. Fujita, K. Koyama, R. Kusakabe, K. Minami, H. Inoue, S. Nishizawa, M. Tsuji, T. Nishiki, M.
Hori, L. Maddegedara, and N. Ueda, “Fast Scalable Implicit Solver with Convergence of Physics-Based
Simulation and Data-Driven Learning,” Research Poster for SC20, 2020.
[12] J. Nathan Kutz, Steven L. Brunton, Bingni W. Brunton and Joshua L., Dynamic Mode Decomposition:
Data-Driven Modeling of Complex Systems, Society for Industrial & Applied Mathematics, 2016.
[13] The American Society of Mechanical Engineers, ASME V&V 10-2006: Guide for verification and validationin computational solid mechanics, 2006.
[14] Supercomputer Fugaku, RIKEN Center for Computational Science [Online].
https://www.r-ccs.riken.jp/en/postk/project[15] Toshio Yoshida, Fujitsu High Performance CPU for the Post-K Computer, IEEE Hot Chips: A Symposium on
High Performance Chips, San Jose, August 21, 2018.
[16] Yuichiro Ajima et al.,“The Tofu Interconnect D,”IEEE Cluster 2018, 2018.
[17] Oakbridge-CX Supercomputer System, Information Technology Center, The University of Tokyo [Online].
https://www.cc.u-tokyo.ac.jp/en/supercomputer/obcx/service/[18] EBE kernel for A64FX CPUs [Online]. To be available.
[19] G. Karypis and V. Kumar, “A fast and high quality multilevel scheme for partitioning irregular graphs,” SIAM
Journal on scientific Computing, vol. 20(1), 1998, pp. 359–392.
[20] G. H. Golub and Q. Ye, “Inexact preconditioned conjugate gradient method with inner-outer iteration,” SIAM
Journal on Scientific Computing, vol. 21(4), pp. 1305–1320, 1999.
[21] Strong ground motion of The Southern Hyogo prefecture earthquake in 1995 observed at Kobe JMA
observatory, Japan Meteorological Agency, [Online].
https://www.data.jma.go.jp/svd/eqev/data/kyoshin/jishin/hyogo_nanbu/dat/H1171931.csv②シミュレーション内で生成されるデータの学習による反復 solver 性能の向上- 69 -

物流・エネルギーの最適化など都市のスマート化が supercomputing のフロンティアの一つと考えられているが (e.g. [1])，これらに加えて今回の新型コロナウイルス感染症（COVID-19）のような都市・社会への脅威に対する安全性・頑健性の検討においても supercomputing の貢献が求められている。例えば、
COVID-19 だけではなく COVID-19 と自然災害のような複合的・致命的な脅威が生じることが懸念されており、従来よりもさらに精緻な都市の状態推定が社会から要請され、これに応えることが喫緊の課題となっている。ひとつの方法として、Big Data & Extreme Computing (BDEC) [2]を前提に、都市の大量の計測データと同化させたリアルタイムで更新される超詳細な都市のデジタルツインを構築し、これを用いた超高詳細なシミュレーションを行うことで、都市の状態推定の高精度・高分解能化を図ることが考えられる。しかしながら、都市は、複雑かつ巨大であるため、この超高詳細なシミュレーションの解析コストは膨大となってしまい、上記を実現するためには解析コストの削減が渇望されている状況にある。
本研究では、都市の脅威となり得る災害の一つである地震を例にこの解析コストの軽減を試みた。な
お、都市の地震問題は、数理的には非線形時間発展問題、計算科学・計算機科学的にはランダムメモリ
アクセス卓越型の低次有限要素解析という性質の悪い難題であり、革新的な開発が必要となる。本研究では、HPC と AI の融合をさらに拡張した新しいアプローチである equation-based modeling と data-driven
learning の HPC による融合により、新たな手法を開発した。具体的には、equation-based modeling の求解過程の解析空間で生じるデータを用いた data-driven learning により局所的な誤差に対する解空間の拘束を与えることで解の収束過程を高速化する、新しい HPC ベースの Fast scalable implicit solver である IRIS を開発し、
上述のように性質が悪い都市地震の問題で、従来法や current state of art である GAMERA [3] (SC14の Gordon Bell Prize finalist) と比較して格段に解析コストが削減されることを示した。これは、従来手法である PCGE と比べて time-to-solution では 14.3 倍、state of art である GAMERA と比べて time-to solution
で 7.52 倍の性能改善である。IRIS を用いることで、図 2.2.1 に示すような従来は難しいとされてきた超高
詳細な都市モデルを用いた解析が可能となった。本研究は、物理空間・サイバー空間・解析空間内の Big
Data が equation-based modeling と data-driven learning による HPC に基づく Extreme Computing によって真価を発揮することで従来の HPC 以上の付加価値を創成し、都市解析のフロンティアを拡げるとともに、supercomputing の都市・社会の守護者としての能力をさらに後押しする innovation の一つになると期待される。
大地震時の都市挙動を信頼性高く評価するためには、超高詳細な都市のデジタルツインの挙動解析を行う必要があり、局所的に超詳細・複雑な構造がある広い領域での媒体の非線形時間発展問題を扱うこととなる。具体的には、解析対象領域は、103×103×101～2 m の範囲であり、かつ、構造物が 10-2～1 m 程
度の非常に複雑な幾何形状をもつ。そのため、応力フリーの境界条件を解析的に満足し、計算コストを抑えつつこのような複雑な形状のモデル化に優れた low order solid element（非構造四面体二次要素）を用いた 3 次元動的非線形有限要素法により連続体の非線形時間発展問題を倍精度で解くこととなり、各時
間ステップ n において An δun=fn なる方程式を 103～4 時間ステップ分解くこととなる。δun の次元は 1010～程度の巨大な問題であり、また、An は疎行列であるものの非構造低次要素であるため演算量と比べてランダムアクセスコストが卓越するという最近の計算機アーキテクチャに親和性が低く高効率な解析を行うことが難しいという性質がある。以上を踏まえて、An δun=fn の求解の収束性を改善するための方法として、current time step 以前の解析結果を用いた data-driven によって高次モードまでの収束過程を学習した
物理シミュレーションのサロゲートモデルを構築し、これを前処理として用いることで収束性を改善す- 70 -

るアルゴリズムを開発した。このアルゴリズムでは、収束性が改善されるとともに、An に関するランダ
ムアクセスの演算が data-driven の連続アクセスの演算へ変換される効果があり、最近の計算機アーキテ
クチャにおいても効率的な演算が可能となる。なお、このアルゴリズムは、一般の収束性の悪い非線形
動的問題に適用可能な汎化性に優れた手法である。本手法によりサロゲートモデル主体の演算に演算量が移行されるが、An に何らかのベクトルをかける演算が主要カーネルとして残るため、これを最新のアーキテクチャにあわせて効率的なチューニングを行うことも本質であり、あわせて最新の CPU アーキテ
クチャにあわせたチューニングも行うことで、前述のように、current state of art と比較して格段の高速化を達成している（図 2.2.2 及び図 2.2.3 参照）。

図 2.2.1 Computation example of 11,002,859,706 degrees of freedom; 2,593,545,132 element model of fully
coupled ground and aboveground/underground structure earthquake shaking analysis using IRIS on 98,304
CPU cores (2,048 nodes) of Fugaku: (a) The 2 × 2 km domain was modeled with 0.5-m 2nd ordered
tetrahedral elements with 9323 buildings and (b) an underground complex with 645 tenants. (c) The complex
nonlinear seismic response was computed as a result of the three-dimensional configuration.

- 71 -

図 2.2.2 Size-up scalability on Fugaku. Numbers in brackets indicate ratio to FP64 peak FLOPS.

図 2.2.3 Speedup scalability on Fugaku.
[1] HPC Connects Plenary: The Century of the City, SC17: International Conference for High Performance
Computing, Networking, Storage and Analysis, 2017.
[2] Big Data and Extreme-scale Computing. [Online].
https://www.exascale.org/bdec/sites/www.exascale.org.bdec/files/whitepapers/bdec_pathways.pdf[3] T. Ichimura, K. Fujita, S. Tanaka, M. Hori, M. Lalith, Y. Shizawa, and H. Kobayashi, ``Physics-based urban
earthquake simulation enhanced by 10.7 BlnDOF x 30 K time-step unstructured FE non-linear seismic wave
simulation,'' in Proceedings of the International Conference on High Performance Computing, Networking,Storage and Analysis (SC'14), IEEE Press, 2014, pp. 15-26.

- 72 -

③ニューラルネットワークによりデータ構造を組み換えた効率的な前処理の高度化③-(a) はじめに
近年 HPC 分野においては、従来からの equation-based method に基づくアプリケーションだけでなく、
data-driven method に基づくものや、equation-based method と data-driven method の組み合わせに基づくものなど、アプリケーションの多様性が増している。また、x86/Arm/Power などの各種 CPU、NVIDIA/AMD/Intelなどの各種 GPU など、commodity から supercomputer systems に至るまで幅広い階層において計算機アー
キテクチャの多様性も増している。これらの多様なアプリケーションにおいて多様な計算機アーキテク
チャを有効利用するため、少ない開発コストで複数の計算機システムにおいてアプリケーションを実行可能とし(portability)十分な実行性能(performance portability)を実現することが課題となっている。
こ の よ う な 背 景 の も と 、 CPU コ ー ド に 指 示 子 を 加 え る こ と で 複 数 の ア ー キ テ ク チ ャ 間 で
portability/performance portability を実現する方法として directive-based parallel programming model が開発されてきた。これまで directive-based parallel programming model の一種である OpenACC [1]を用いて多数の
equation-based アプリケーションを GPU にポーティングした例が報告されており(e.g., [2-6])、native なプ
ログラミング環境(i.e., CUDA for NVIDIA GPU)を使った場合と実用上遜色ない性能が実現できている例も報告されている(e.g., [7])。その一方で、今日開発が加速している data-driven method と equation-based
method を組み合わせたより多様なアプリケーションに対して directive-based programming model の有効性が示された例はまだ少ない。そこで本研究では data-driven method と equation based simulation を融合した
アプリケーションを OpenACC により GPU に porting することで、近年の多様化するアプリケーションにおいても directive based parallel programming が有効であることを示す。
本研究では、data-driven method を equation-based method に組み込んで partial differential equation (PDE)
solver を高速化する、「富岳」[9]用に開発した neural network (NN)-accelerated implicit solver [8]を対象とする。通常 NN などの data-driven method を equation-based method と組み合わせる際には NN の推定精度不足が問題となるが、この手法では対象となる PDE の特性を反映するグリーン関数 (GF)を介して NN を活用する前処理方法(GF-based NN preconditioner)を開発することで、data-driven method の精度の問題を解決している。PDE solver に含まれる疎な計算・ランダムデータアクセスが、近年の計算機と親和性が高い密な計算・連続データアクセスとなる NN に移ることで、「富岳」等の CPU ベースのシステムにおいて従来の PDE ベース solver 比で 4.26 倍高速化している。また、これらの計算は空間的にローカルとなるため、「富岳」の高並列環境上で高い scalability も実現している。この solver algorithm は近年の幅広い計算機アーキテクチャにおいて有効であると考えられ、GPU アーキテクチャの特性を踏まえたカーネルアルゴリ
ズム・実装を構築することで GPU においても高効率計算が実現できると期待される。
以下の構成を示す。Section ③-(b)において対象となる solver を説明し、Section ③-(c)において solverの GPU へのポーティング方法を説明する。ここでは、スケーラブルな solver algorithm をベースに GPU
アーキテクチャに併せてカーネルアルゴリズムを構築し、OpenACC を用いてアプリケーションを実装する。Section ③-(d)において、開発手法の性能を計測する。ここでは、CPU 用のカーネルアルゴリズムをそのまま GPU で実行した場合と GPU に適したカーネルアルゴリズムを構築した場合の性能差や、高実装コストだが自由度の高い CUDA を活用した実装との性能差を比べることで、OpenACC のような低実装- 73 -

コストな実装方法を使った場合でも GPU アクセラレータにおいて高い性能を実現できることを示す。
Section ③-(e)でまとめを示す。
③-(b) Solver with Green’s function-based NN preconditioner
本研究では、neural network (NN)-accelerated implicit solver の一例として、対象とする PDE の特性を反映した Green's function を介して NN を構築し反復解法の前処理として用いる Green's function-based NN
preconditioned implicit solver [8]を対象とする。以下、target とする PDE, GF-based NN の構築方法、及び、
この NN の implicit solver 中での活用方法を説明する。
③-(b) A. ターゲット問題
本研究では、PDE ベースの物理シミュレーションの一例として、非均質物性領域中の波動伝播を対象
とする。
𝜕2 𝑢

𝜕

𝜕𝑢

𝜌 𝜕𝑡 2𝑖 = 𝜕𝑥 (𝑐𝑖𝑗𝑘𝑙 𝜕𝑥𝑘 ) + 𝑓𝑖 . (2.3.1)𝑗𝑙

ここで、𝑥、𝑡、𝜌、𝑢、𝑐、𝑓は位置、時間、密度、変位、剛性テンソル、外力を示す。式(2.3.1)の離散化に
1

1

4

2

は簡単のためボクセル有限要素法を用い、Newmark-𝛽法 (𝛽 = , 𝛾 = )による implicit time integration を用いる。ターゲット問題は、
𝐴𝛿𝑢 = 𝑓, (2.3.2)となる。ここで、
𝐴=

4
𝑑𝑡 2

𝑀+

2
𝑑𝑡

𝐶 𝑛 + 𝐾 𝑛 , (2.3.3)

𝑓 = 𝑏𝑛 − 𝑞𝑛−1 + 𝐶 𝑛 𝑣 𝑛−1 𝐴 (𝑎𝑛−1 +4𝑑𝑡𝑣 𝑛−1 ). (2.3.4)

であり、𝑀、𝐶 𝑛 、𝐾 𝑛 、𝛿𝑢は質量マトリクス、減衰マトリクス、剛性マトリクス、変位増分である。上添
え字の𝑛は時間ステップを示す。式(2.3.2)を解くことで𝛿𝑢を求め、
𝑢𝑛 = 𝑢𝑛−1 + 𝛿𝑢,𝑣 𝑛 = −𝑣 𝑛−1 +𝑎𝑛 = −𝑎𝑛−1 −

2
𝛿𝑢,
𝑑𝑡

4 𝑛−1
4
𝑣
+ 2 𝛿𝑢,𝑑𝑡
𝑑𝑡

𝑞𝑛 = 𝑞𝑛−1 + 𝐾 𝑛 𝛿𝑢.
により変位・速度・加速度・内力を更新することとなる。本研究で対象とする問題は 109 自由度以上の大規模問題となる。そこで、未知数が 109 程度となる式(2.3.2)の求解を高速化することがターゲット問題となる。
③-(b) B. GF-based NN predictor
NN 等の data-driven method は必ずしも結果の精度保証がされていないため、本研究では NN を式(2.3.2)の直接求解には使わず、反復解法の前処理における線形方程式の解の predictor として用いる。この solver
が通常の equation-based method に基づく solver に対して高速となるためには equation-based method におけ- 74 -

る前処理よりも効率的に精度を稼ぐ(i.e., 精度あたりのコストを下げる)ことが重要となる。NN による
predictor として一番 straight forward な方法は、系全体の挙動を追うことが可能なサロゲートモデルを NNを用いて構築する方法であり、このアプローチは小規模問題において有効であることが示されている[10][11]。その一方で、大規模問題においてはモードの数が増えるだけでなく、各モードを表わすデータ(i.e., (𝑥, 𝐴𝑥)の組)も大規模になるため、系全体の挙動を追おうとするアプローチで推定性能をあげるためには多大なコストがかかってしまうこととなる。そこで、PDE の特性を反映するグリーン関数（GF）を介してローカルな NN を構築することで低コストに高精度な推定解を得る。以下、GF-based NN の構築方
法の概略を説明する。
PDE
𝐿(𝑥)(𝑎(𝑥)) = 𝑏(𝑥), (2.3.5)
に対し、GF 𝑔(𝑥, 𝑠)は以下の式を満たす。
𝐿(𝑔(𝑥, 𝑠)) = 𝛿(𝑠 − 𝑥), (2.3.6)
ここで、𝐿、𝑏、𝑎、𝛿 は linear differential operator, 既知の分布、システム応答、及び、Dirac delta functionを示す。この GF を用いることで、システムの応答は
𝑎(𝑥) = ∫ 𝑔(𝑥, 𝑠)𝑏(𝑠)𝑑𝑠. (2.3.7)と求めることができる。このように GF は PDE の特性を反映するもので、畳み込み計算により線形方程式を直ちに求解できる。その一方で、離散化された GF の分布を求める計算は、PDE を離散化した線形方程式を解く計算と同等のコストがかかるため、GF を使ったアプローチは通常使われない。そこで本研究では、以下のように GF を NN にて高速に推定することで GF のコストの問題・NN を利用する際に課題となる精度の問題を解決する。
ターゲット問題は非均質物性領域であるため、GF の分布も物性分布に応じて変化することとなる。離
散化した GF は図 2.3.1 にあるように複雑な分布をしており、少量のデータによって高次モードまでを効率良く推定することが難しい問題となっている。そこで、複雑形状を持つ非均質領域の GF が、均質領域の GF の分布の簡単な変換を用いて近似できると仮定し、以下の 2 段階を経て非均質領域の GF を推定する。
まず、非均質物性を代表する物性値を一つ選び、この物性の均質領域における離散化された GF を求める。波動方程式の GF は距離の二乗に反比例して減少するため、GF が有意な値を持つ𝑁𝑥 × 𝑁𝑦 × 𝑁𝑧 節点の領域での GF 値を求める（それより外側では GF の値は 0 と近似する）。
ここでは、ランダムな入力場𝑥と、それに𝐴をかけた𝐴𝑥の組を使い、数値最適化により𝑁𝑥 × 𝑁𝑦 × 𝑁𝑧 節点の GF の分布を求める。
次に、均質な物性領域における GF を不均質物性領域における GF にマッピングする NN を構築する。
ここでは、GF の 9 成分𝐺𝑖𝑗𝑁𝑁 (𝑖, 𝑗 = 1,2,3)が以下の 8 つの係数𝑐𝑖 (𝑖 = 1,2, … ,8)により展開できると仮定し、𝐺𝑖𝑖𝑁𝑁 = 𝐺𝑖𝑖𝑏𝑎𝑠𝑒 (𝑐1 + 𝑐2 𝑥 + 𝑐3 𝑦 + 𝑐4 𝑧), (2.3.8)
𝐺𝑖𝑗𝑁𝑁 = 𝐺𝑖𝑗𝑏𝑎𝑠𝑒 (𝑐5 + 𝑐6 𝑥 + 𝑐7 𝑦 + 𝑐8 𝑧) for 𝑖 ≠ 𝑗 , (2.3.9)
これらの値を推定する NN(input: 𝑁𝑥−1 × 𝑁𝑦−1 × 𝑁𝑧−1要素の物性値、output: 𝑐𝑖 (𝑖 = 1,2, … ,8))を構築する。
ここでは、fully connected 8-4-4-4-4-4-4-25 feed-forward network を用いた Classifier NN を用いる(実装にはChainer [12]を用いた)。
これにより、非均質物性領域の GF を直接推定する場合と比べて、少ないデータ量・パラメータ数で高精度な GF を推定する NN を構築可能となっている。図 2.3.1 に学習に 16.2 million 組のデータを使って構- 75 -

築した NN を用いて推定した GF の分布を示す(𝑁𝑥 = 𝑁𝑦 = 𝑁𝑧 = 7; input の物性パラメータ数は 63=216、GF の出力値は 73×9 成分となる。
図の上段が正確な GF で、下段が NN にて推定された GF の誤差である。
ここでは、Section③-(d)における性能計測問題の離散化設定における GF をプロットしている。図から、
GF の各成分がそれぞれ複雑な高次モードを有しているにも関わらず、本手法により複数の物性を持つ領域における GF の各成分を高い精度で求められていることがわかる。
求めた𝐺𝑖𝑗𝑁𝑁 (𝑖, 𝑗 = 1,2,3)により𝐴𝑧 = 𝑟の解𝑧は
𝑁

𝑁

𝑁

𝑁

𝑁

𝑁

𝑧𝑖 (𝑖𝑥 , 𝑖𝑦 , 𝑖𝑧 ) = ∑3𝑗=1 ∑𝑗𝑧𝑧=1 ∑𝑗 𝑦=1 ∑𝑗𝑥𝑥=1 {𝐺𝑖𝑗𝑁𝑁 (𝑗𝑥 , 𝑗𝑦 , 𝑗𝑧 )𝑟𝑗 (𝑖𝑥 + 𝑗𝑥 − 2𝑥 , 𝑖𝑦 + 𝑗𝑦 − 2𝑦 , 𝑖𝑧 + 𝑗𝑧 − 2𝑧 )} for 𝑖 = 1,2,3
𝑦(2.3.10)と近似する。ここで、𝑧𝑖 (𝑖𝑥 , 𝑖𝑦 , 𝑖𝑧 ) はボクセル節点𝑖𝑥 , 𝑖𝑦 , 𝑖𝑧 における第𝑖成分の値、𝐺𝑖𝑗𝑁𝑁 (𝑗𝑥 , 𝑗𝑦 , 𝑗𝑧 )は𝑗𝑥 , 𝑗𝑦 , 𝑗𝑧 における GF の𝑖, 𝑗成分の値を示す。

図 2.3.1 Accuracy of GFs obtained by NNs. Upper row indicates the true values, while the lower row indicatesthe errors of the NN-approximated values.
③-(b) C. Scalable solver algorithm using GF-based NN predictor
可変前処理付き共役勾配法[13]に GF-based NN を組み込むことで scalable な solver algorithm を構築する。
通常の前処理 CG 法では𝐴−1を近似する固定マトリクス𝑀を用いて探索方向𝑧 = 𝑀𝑟を計算するが、可変前処理付き CG 法では𝐴𝑧 = 𝑟を近似的に解いた𝑧を探索方向に用いる。これにより、inexact な計算となる
GF-based NN predictor を前処理内で使うことが可能となる。以下、Algorithm 1 に沿って説明する。

- 76 -

Preconditioner: GF の形状は内部領域と領域境界では大きく変わるため、境界から十分離れた領域ではGF-based NN predictor を用いて𝑧 = 𝐴−1 𝑟を推定し(Algorithm 1 line 8)、境界付近は通常の CG 法 (3×3 block
Jacobi preconditioned CG 法)で粗く解く(Algorithm 1 line 9)。なお、MPI 領域分割における領域境界においては GF 計算に必要な通信を回避するため、境界付近と同様に CG 法でアップデートする。これらの境界
領域のアップデート時には、省メモリでスケーラビリティに優れた 3×3 block Jacobi Preconditioned CG 法を用いる（図 2.3.2）。
Matrix-free matrix-vector product: 限られたメモリ量を有す計算機において大規模計算を実施するため、メモリ使用量の多い全体行列𝐴はメモリに格納せずに、行列ベクトル積𝑓 = 𝐴𝑢を計算する度に要素マトリ
クス𝐴𝑒 と要素右辺ベクトル𝑢𝑒 との積を計算し、全体左辺ベクトル𝑓に足しこむ Element-by-Element 法[14]を用いる。全体行列を格納する場合と比べてメモリ使用量だけでなくメモリ転送量も削減可能なため、メモリバンド幅に対し演算性能が相対的に高い計算機向けのアルゴリズムになっている。Algorithm 1 line

- 77 -

16 における行列ベクトル積とともに、Algorithm 1 line 9 の境界部分のアップデートの CG 法内での行列ベクトル積においても EBE 法を用いる。
Mixed precision arithmetic: 計算の最終的な出力は倍精度が必要であるものの、前処理の計算は inexact
で良いため、前処理計算(Algorithm 1 line 7-10)は単精度で実施する（前処理以外の計算は全て倍精度で実施）。高いピーク性能のハードウェアが使えるだけでなく、メモリへの読み書きや隣接通信時のデータ移
動量も半分になる。

図 2.3.2 Using GF-based NN predictor. White parts are predicted by GF-based NN while the colored parts
near inter-process boundary and physical domain boundaries are solved roughly using a conjugate gradientsolver. nef=3 is used in this study.
③-(c) GPU porting of solver with Green’s function-based NN preconditioner using OpenACC
GF-based NN preconditioned solver algorithm は NN による密な計算＋構造格子となるため、GPU を含む
幅広いアーキテクチャにおいて高い性能が期待される。その一方で、カーネルレベルでみると GPU において性能を出しやすいアルゴリズムとそうでないアルゴリズムがあるため、カーネルレベルにおいて適切なアルゴリズム選択が必要となる。そこで本研究では GPU メモリからの読み書きの削減や coalesce 化による帯域の有効活用、ローカルなカーネル内で GPU アーキテクチャに適したカーネルアルゴリズムを構築した。また、GPU においてはその高い計算能力のために通信がボトルネックに陥りやすいという傾向にあるため、MPI process 間での通信削減も実施する。なお、前提として、全ての計算を GPU にオフロードし、CPU は GPU の管理のみに利用するよう OpenACC 化している。以下、主要カーネルとなる、- 78 -

GF-based NN predictor、Matrix-vector product kernel、Communication kernels についての GPU に適したアルゴリズム構築・カーネル実装を説明する。
Avoiding data recurrence in matrix-vector product kernel: 行列ベクトル積カーネルの並列計算においては左辺ベクトルへの足し込み時において生じる data recurrence を回避する必要がある。一番単純な方法は図2.3.3 の左図のように要素を 8 つの色にカラーリングし、それぞれのカラー内での data recurrence を回避する方法である。この方法はシンプルで多数のスレッドを実行する GPU においても活用できるが、stride 2のデータアクセスが生じることとなる。近年の GPU には hardware accelerated atomics が実装されており、特に今回対象とする A100 GPU といった NVIDIA の GPU ではグローバルメモリで高い atomic 演算のスループットを実現しているため、これを活用した data recurrence の回避により性能向上が期待できる。図2.3.3 右図にあるように、atomic add を活用することで、カラーリングアルゴリズムでの stride 2 のデータアクセスが連続アクセスに置き換えられており、高速な atomic 演算が可能な計算機においては性能向上が期待される。

図 2.3.3 Avoiding data recurrence in the matrix-vector product kernel. a) Left: Using coloring of elements, b)
Right: using atomic add.

- 79 -

Element wise computation in matrix-vector product kernel:近年の GPU はスレッドあたり大量のレジスタを有しており、レジスタ上での計算量を増やしてでも data access を減らした方が性能が出しやすい傾向となる。そこで、レジスタ数が相対的に少ない CPU 等を対象とした一般的な実装ではあらかじめ計算しメ
モリ上に格納しておく中間変数においても、GPU 上ではオンザフライで都度生成するようにカーネルア
ルゴリズムを変更する。
これにより、
行列ベクトル積カーネル内での GPU メモリ読み出し量が削減され、速度向上が期待される。
GF-based NN predictor kernel: 本カーネルは、周辺グリッドの情報をローカルな節点に畳み込む計算となり、図 2.3.4 のような 3 重ループ内に 3 重ループがあるような単純な構造となっている。GPU における多数のスレッドを有効活用するため、外側の 3 重ループを collapse することで並列化にかかるループ長を十分長く取るとともに、GPU で高いデータアクセス性能を実現するためカーネル内で読み出す変数の配列次元を適切に設定することでデータアクセスは全て coalesced access とする。また、外側ループをアンロールすることで、最内ループでロードする変数を再利用し broadcast load 数を減らすことでさらなる性能向上も可能である。

図 2.3.4 GF-based NN predictor kernel- 80 -

Mapping of processes for efficient communication: 現代の GPU 搭載計算機は複数の GPU と複数の通信ポートを搭載していることが多く、これらの構成にあわせたプロセス配置や通信ポート割り当てによりノード間通信を減らし、通信バンド幅を有効活用するために重要となる。ここでは、ノード内の GPU 間の通信がノード間の通信よりも高速で、かつ、ノード内に搭載された GPU と通信ポート間の転送速度に偏りがある環境を考慮し、図 2.3.5 のようなプロセスマッピングを採用する。この際、各 GPU に近い通信ポートを使いつつ、各通信ポートにおける通信量のバランスが保たれるようプロセス配置している。上記は本研究での性能計測に用いる 2 次元プロセス分割・AI Bridging Cloud Infrastructure （ABCI）の計算
ノードの構成を例に作成したマッピングだが、他の構成や 3 次元プロセス分割においても同様に適切なマッピングが可能である。

図 2.3.5 Mapping of processes for efficient communication. Here, mapping of processes for a systemcomprising compute nodes each with 8 GPUs and 4 NICs are shown. By using b) instead of a), maximum
inter-node communication volume can be reduced by 40% per NIC.
③-(d) 性能計測③-(d) A. Problem used for measurement
非均質物性領域における大規模波動伝播解析を想定し、性能検証の継続性が重要であることから、[8]での性能計測問題と同じく、CT スキャンにより作成した頭部モデル[15]における波動伝播解析の性能を計測する(図 2.3.6)。この計測の問題設定は本開発手法の地震問題での有効性を十分検証できるものであり、実際に本開発手法をベースとした手法を超高詳細断層-都市地震解析（HPCAsia2022）に適用している。
ここでは、𝑉𝑠 =50～120 m/s と可変とし、他の物性値は damping ℎ=0.001、density 𝜌=1,000kg/m3、𝑉𝑝 =200m/s
で固定とする。𝑑𝑡=0.001s で離散化し、512×512×512 elements (513×513×513 節点、節点あたり x,y,z 方向の- 81 -

3 自由度の合計 405,017,091 自由度)。Weak scaling 計測時はこれを x,y 方向に敷き詰める。本来は時系列問
題を解くこととなるが、時間方向に行列の性質は変わらず反復法 solver の収束特性は大きく変わらない
ため、本研究では入力としてホワイトノイズを𝑓に設定し、底辺ディリクレ固定に対する 1 タイムステッ
プの求解にかかる時間を計測する。NN の構築コストは小さく、あらかじめ構築しておくことができるた
め、NN の推論を含むソルバー部のみを計測対象とする。

図 2.3.6 Head model used for performance measurement. Young's modulus and von Mises stress response of
the model under gravity is shown.
本研究では、
比較対象の solver として超並列環境において大規模な波動問題を解く際に広く使われる、
equation based solver の一種である CGBJ (3×3 ブロックヤコビ前処理付き共役勾配法)を用いる。CGBJ は
節点毎にローカルに求めた逆行列を前処理に用いているためロードバランスがとりやすく、また計算が行列ベクトル積のみからなっているため、行列ベクトル積をオンザフライで実施する Element-by-Element法と組み合わせることで省メモリで実装可能であり、超大規模問題の求解に適した方法となっている。
本計測で用いる CGBJ solver は提案手法で開発した solver における前処理部(Algorithm 1 line 7-10)を 3×3
ブロックヤコビ法(i.e., 𝑧 = 𝐵𝐽−1 𝑟 により探索ベクトル𝑧を設定。ここで、𝐵𝐽 : 3×3 ブロックヤコビ行列)に代えることで実装しているため、提案手法と同じチューニングレベルでの比較対象となっている。
③-(d) B. Performance measurement environment
近年多数のスーパーコンピュータシステムで用いられている NVIDIA A100 [17]及び V100 [18] GPU を搭載した産総研の ABCI[16]を用いて性能計測を実施する。表 2.3.1 に計算ノードの構成を示す。システ- 82 -

ムは A100 を搭載する計算ノード(A)120 台、V100 を搭載する計算ノード(V)1,088 台からなり、計算ノー
ド(A)は全ノードフルバイセクションバンド幅、計算ノード(V)は 34 ノードからなるラック内においてフルバイセクションバンド幅で相互接続される。計算ノード(A)においては CPU と GPU の FP64 ピーク性
能比は 14.0 倍（メモリバンド幅は 30.4 倍）、計算ノード(V)においては CPU と GPU の FP64 ピーク性能
比は 10.2 倍（メモリバンド幅は 14.1 倍）となっており、A100 ノードにおいては GPU のメモリバンド
幅が特に増強されていることがわかる。GPU 計測においては nvhpc 21.2、cuda 11.2.2、openmpi 4.0.5 を用い、コンパイラオプションとして-fopenmp -fastsse -O3 -Minline=levels:10 -Mcuda=cuda"11.2",ptxinfo -acc
-ta=tesla:cc"80",loadcache:L1,fastmath -Minfo=accel -mcmodel=medium -Mlarge_arrays、実行時オプションとして--mca btl_openib_want_cuda_gdr 1 -x UCX_MEMTYPE_CACHE=n -x UCX_MAX_EAGER_RAILS=1 -x
UCX_MAX_RNDV_RAILS=1 -x PGI_ACC_BUFFERSIZE=50M を 用 い る (CUDA-Aware MPI and GPU
Direct RDMA [19] communication を用いる)。CPU 計測においては intel 2020.4.304、intel-mpi 2019.9 を用い、コンパイラオプションとして-O3 -xCORE-AVX512 -qopenmp -qopt-report -qopt-zmm-usage=high を用い、OpenMP においてノード内の全コアを利用する(計算ノード(A)では 9 コア/プロセス、計算ノード(V)では 10 コア/プロセス)。以下の計測では、1 プロセス (1 GPU)あたり 256×256×512 要素に問題サイズを固定して計測する。
表 2.3.1 Configuration of ABCI system③-(d) C. Solver performance on GPU-based system
本手法の GPU 上での実行性能を確認する。まず、カーネルレベルでの性能を計測する(表 2.3.2)。
EBE kernel: 本カーネルは中間変数が多くデータアクセスも多いため、レジスタ数が限られている CPU
においては spill/fill が発生し、ピーク性能は 11.2%であったが、GPU においては多数のレジスタにより
spill/fill が回避され、データアクセスコストが削減されるため、高いピーク性能比が実現した。まず、図2.3.3a に示すカラーリングの場合、EBE kernel の性能は 28.2 ms、5.33 TFLOPS (FP64 ピーク性能比の 55%)となり CPU と比べて高い性能となっているが、図 2.3.3b の atomics を有効利用するアルゴリズムを用いた場合の性能は 21.5 ms、6.99 TFLOPS (FP64 ピーク性能比の 72%)と、カラーリング時のさらに 1.31 倍速となった。さらに、中間変数をメモリから読みださずに都度再計算するアルゴリズムを利用することで、- 83 -

CPU 比で 93.7 倍の高速化(ピーク性能比 75.0%)につながっている。なお、OpenACC において routine ディ
レクティブを用いることによって parallel 領域内部から他のデバイス関数を呼び出すことが可能であるがオーバーヘッドが大きくスレッドの細粒度の制限があるため関数呼び出しは手動で inline 展開している。
CPU 実装では、メインのループ内部に中間変数計算のための行列-行列積関数などの関数呼び出しが複数含まれており、更に一部関数ではその内部で更に関数呼び出しを行うような状態となっていた。OpenACC
でも routine ディレクティブを用いることによって parallel 領域内部からも他のデバイス関数を呼び出すことが可能となっているが、呼び出しのオーバーヘッドに加えてスレッドの細粒度の制限などにより、性能を出し切ることが難しくなる。
そこで今回は parallel 領域内部の関数呼び出し全てに対して手動で inline
展開を行う。inline 展開によってレジスタ消費量は増加するが、今回の実装ではレジスタスピルは確認されず、更なる速度向上が期待できる。
GF-based NN predictor kernel: もともと演算密度が高くデータアクセスの少ない区間であるため、CPU
では FP64 ピーク性能比の 70.5% (FP32 ピーク性能比の 35.2%)と高い性能が出ているが、A100 GPU においてはさらに高い性能(FP32 ピークの 49.4%, FP64 ピークの 98.8%に相当する 9.63 TFLOPS)がでていることがわかる。これも多数のスレッド・レジスタによりデータアクセスや演算レイテンシが適切に隠蔽されているためと考えられる。
このように、CPU 実装と比べると、EBE kernel、Predictor kernel ともに大幅に高速化していることがわかる。これらのカーネルレベルの高速化により、ソルバー全体で A100 ノードで CPU 比 38.9 倍速、 V100
ノードで CPU 比 22.4 倍速となっている(表 2.3.2)。
実装コストが高いがプログラミングの自由度が高い CUDA における実装との性能比を比べると、A100
ノードにおいては Outer EBE kernel と Predictor kernel でそれぞれ対 CUDA 比で 93.0%、85.5%の性能であった。また、V100 ノードにおいては Outer EBE kernel と Predictor kernel でそれぞれ対 CUDA 比で 91.8%、
82.7%の性能であった。このように GPU に適したアルゴリズムを設計することで OpenACC による実装がCUDA による実装と比べて実用上遜色ない性能となっていることがわかる。
表 2.3.2 Performance of proposed solver on compute node (A)× 1 node

- 84 -

表 2.3.3 Performance of proposed solver on compute node (V)× 2 nodes

図 2.3.7 において GF-based solver の性能を示す。CGBJ は節点毎にローカルに求めた逆行列を前処理に用いているためロードバランスがとりやすく、また計算が行列ベクトル積のみからなっているため、行
列ベクトル積をオンザフライで実施する Element-by-Element 法と組み合わせることで省メモリで実装可能であり、超大規模問題の求解に適した方法となっている。図から、CGBJ が求解に 42 反復要しているのに対し、本手法では 5 反復で相対誤差|𝐴𝛿𝑢 − 𝑏|/|𝑏| < 10−8まで解けていることがわかる。ここから、本手法の前処理の精度が従来からの equation-based な前処理と比べても高いことがわかる。
なお、
GF-based
NN は初期解推定においても使うことが可能であり、実際、図 2.3.7 の GF-based NN initial solution + CGBJのように初期解の精度を改良することが可能である。ただし、解の refinement において GF-based NN を使わない場合、性能改善の幅は小さく、本手法のように前処理内において解の refinement を実施する際にGF-based NN を使うことが性能向上に重要となっていることがわかる。本手法・CGBJ ともに SIMD
intrinsics により AVX-512 にてチューニングを施した CPU 実装を用いて本問題を CPU（Intel Xeon Gold
Platinum CPU (Cascade Lake))上で計測した場合、GF-based solver の 1 反復あたりの所要時間は CGBJ 1 反復の約 1.64 倍に抑えられ、solver 全体で 4.26 倍の高速化につながっている[8]。高精度な predictor により反復数が大幅に減っていること、
及び、
ピーク性能が出しやすい NN によりピーク性能比が CGBJ の 10.9%から提案手法において 24.1%に向上していることで一反復あたりの実行時間が短く抑えられていることが大幅な高速化につながっている。

- 85 -

図 2.3.7 Comparison of the performance of solver algorithms. Elapsed time on compute node (A) system isshown in the horizontal axis.
また、提案手法は従来手法である CGBJ を GPU で実行した場合と比べて、A100 ノードで 2.54 倍速、V100 ノードで 2.68 倍速となっている。ここで CGBJ の実装は提案手法の EBE をそのまま用いたものとなっており、ピーク性能は A100 ノードで FP64 ピーク性能比の 41.9%とかなり速い実装に相当するが、それでもなおこの高速化が達成されていることから提案手法が有効なアルゴリズムとなっていることがわかる。
③-(d) D. Weak scaling on GPU-based system
最後に、提案手法の大規模問題でのスケーラビリティを測る。ここでは、1 プロセス (1GPU)あたり
256×256×512 要素に問題サイズを固定し、プロセス数と問題規模を拡大する。表 2.3.4 に Compute Node (A)1 ノード(8 GPU) to 32 ノード(256 GPU)の性能を示す。表から、point-to-point communication にかかる時間
が 1～4 ノード間で増えているものの、4 ノード以上ではほぼ一定となっていることがわかる。これは、full bisection network の特性で隣接通信が効率よく行われているためと考えられる。結果として、1 ノード(8 GPU)から 32 ノード (256 GPU)までは 83.4%の高い weak scalability となっていることがわかる。多数ノードにおいてもロードバランスが保たれ、精度混合演算による通信量削減型の前処理が効いているために多数の GPU からなる並列計算機上でも高い性能となっていると考えられる。
表 2.3.4 Weak scaling of proposed solver on compute node (A) system.

- 86 -

③-(e) まとめ
本研究では、「富岳」用に開発した equation-based method と data-driven method を融合した HPC アプリ
ケーションの多システムへの展開の一例として、GPU 上で高効率に実行できる Green's function-based
Neural Network solver を directive-based なポーティング（OpenACC）により実現した。ここでは、Green's
function を介して NN を用いることで高精度・低コストで反復解法の前処理を実施するスケーラブルな
solver algorithm において、計算機アーキテクチャに併せたカーネルアルゴリズムを構築し OpenACC を用
いてアプリケーションを実装することで、GPU において高い性能を低いコード開発コストにて実現した。
開発手法は FP64 ピーク性能比の 64.4%という高い性能が得られ、NN を使わない従来からの solver と比べて本手法を用いることで 2.54 倍速の time-to-solution の高速化が得られた。また、Target native な環境である CUDA を使って同じアルゴリズムを実装した場合のコストを OpenACC と比較し、カーネルレベル
で 85%～93%という高い性能が得られることも示した。このように、適切なアルゴリズム構築によりOpenACC のような低実装コストな実装方法を使った場合でもアクセラレータにおいて高い性能を実現できることを示した。本研究で示した方法は OpenACC の他の directive based programming models においても有効であると考えられ、現在様々な分野で開発が進められている data-driven method と equation-based
method を融合した HPC アプリケーションの accelerator-based systems における低コストで開発・実行に貢献すると期待される。

[1] OpenACC, http://www.openacc.org/, last accessed 2021/08/24[2] Chuang, P.-Y., Foertter, F.S.: An Example of Porting PETSc Applications to Heterogeneous Platforms with
OpenACC. In: Accelerator Programming Using Directives. pp. 3–19. Springer International Publishing (2018).
https://doi.org/10.1007/978-3-319-74896-2_1.
[3] Shan, H., Zhao, Z., Wagner, M.: Accelerating the Performance of Modal Aerosol Module of E3SM Using
OpenACC. In: Accelerator Programming Using Directives. pp. 47–65. Springer International Publishing (2020).
https://doi.org/10.1007/978-3-030-49943-3_3.
[4] Liang, J., Hua, R., Zhang, H., Zhu, W., Fu, Y.: Accelerated molecular dynamics simulation of Silicon Crystalson

TaihuLight

using

OpenACC.

ParallelComputing.

99,

102667(2020).

https://doi.org/10.1016/j.parco.2020.102667.
[5] Xue, W., Roy, C.J.: Multi‐GPU performance optimization of a computational fluid dynamics code using
OpenACC. Concurrency Computat Pract Exper. 33, (2020).
[6] Londhe, A., Rastogi, R., Srivastava, A., Khonde, K., Sirasala, K.M., Kharche, K.: Adaptively accelerating
FWM2DA seismic modelling program on multi-core CPU and GPU architectures. Computers & Geosciences.
146, 104637 (2021). https://doi.org/10.1016/j.cageo.2020.104637.
[7] Yamaguchi, T., Fujita, K., Ichimura, T., Naruse, A., Lalith, M., Hori, M.: GPU Implementation of a
Sophisticated Implicit Low-Order Finite Element Solver with FP21-32-64 Computation Using OpenACC. In:
Accelerator

Programming

Using

Directives.

pp.

3–24.

SpringerInternational

Publishing(2020).

https://doi.org/10.1007/978-3-030-49943-3_1.
[8] Ichimura, T., Fujita, K., Hori, M., Maddegedara, L., Ueda, N., Kikuchi, Y.: A Fast Scalable Iterative Implicit
Solver with Green’s function-based Neural Networks. In: 2020 IEEE/ACM 11th Workshop on Latest Advances- 87 -

in Scalable Algorithms for Large-Scale Systems (ScalA). IEEE (2020).
https://doi.org/10.1109/scala51936.2020.00013.
[9] About Fugaku, https://www.r-ccs.riken.jp/en/fugaku/about/, last accessed 2021/08/24[10] Sappl J., Seiler L., Harders M., Rauch W.: Deep Learning of Preconditioners for Conjugate Gradient Solvers in
Urban Water Related Problems. arXiv preprint arViv:1906.06925 [cs.LG] (2019).
[11] Gotz, M., Anzt, H.: Machine Learning-Aided Numerical Linear Algebra: Convolutional Neural Networks for
the Efficient Preconditioner Generation. In: 2018 IEEE/ACM 9th Workshop on Latest Advances in Scalable
Algorithms for Large-Scale Systems (scalA). IEEE (2018). https://doi.org/10.1109/scala.2018.00010.
[12] Chainer, https://chainer.org/, last accessed 2021/08/24[13] G. H. Golub and Q. Ye, Inexact preconditioned conjugate gradient method withinner-outer iteration, SIAM
Journal on Scientific Computing, 21, pp.1305–1320, 1999.
[14] Winget, J.M., Hughes, T.J.R.: Solution algorithms for nonlinear transient heat conduction analysis employing
element-by-element iterative strategies. Computer Methods in Applied Mechanics and Engineering. 52, 711–815(1985). https://doi.org/10.1016/0045-7825(85)90015-5[15] Arayeshnia, A., Keshtkar, A., Amiri, S.: Realistic human head voxel model for brain microwave imaging. In:2017 Iranian Conference on Electrical Engineering (ICEE). IEEE (2017).
https://doi.org/10.1109/iraniancee.2017.7985315.
[16] AI Bridging Cloud Infrastructure (ABCI), National Institute of Advanced Industrial Science and Technology(AIST) [Online]. https://abci.ai/en/about_abci/, last accessed 2021/08/24[17] NVIDIA A100 TENSOR CORE GPU, https://www.nvidia.com/en-us/datacenter/a100/, last accessed
2021/08/24[18 NVIDIA V100 TENSOR CORE GPU, https://www.nvidia.com/en-us/datacenter/v100/, last accessed
2021/08/24[19] NVIDIA GPUDirect, https://developer.nvidia.com/gpudirect, last accessed 2021/08/24④開発された大規模数値解析コード群の実展開を踏まえた解析手法の安定化及び有効性の確認④-(a) 液状化解析手法の安定化④-(a) A. はじめに計算機環境、計算手法の発達に伴って、physics-based の大規模なシミュレーションが可能になってきた。さらに、膨大なデータが収集可能になってきたことに伴って、big data and extreme-scale computing(BDEC) への期待が高まっている。BDEC では、膨大な観測データやそれを学習データにした機械学習と計算科学のシミュレーションが統合され、信頼性の高いシミュレーションが可能になると考えられている。都市のデジタルツイン構想はその例の１つで、BDEC 用の計算システム (BDEC system) 上に高詳細な都市モデルであるデジタルツインを構築し、時々刻々と得られる観測データをもとにモデルの状態を更新しながら、デジタルツイン上で数値解析による仮想実験を行い、その結果を政策などの意思決定に用いる。このような検討では多数回の解析を行うこととなるため、実際に使えるようにするためには大規模な解析であっても安定的に解析が出来ることが必須である（デジタルツインに限らず、解析が不安定化することは問題となるため、安定的に解析できることは本質的に重要である）。そのため、開発され- 88 -

た大規模数値解析コード群の実展開を目指して、解析手法の安定化及び有効性の確認を本研究では行った。
具体的には、本研究では、解析が不安定化しやすく、かつ、大規模な解析が必要とされる液状化シミ
ュレーションを例とすることとし、季節変動によって変化する地下水位にともなって変動する地震リス
クを計算するための手法を開発した。従来の解析では、地下水位は地表面から 2m というように恣意的に設定されることが多かった。本研究では、より正確な地震リスクを行うことを目指し、浸透流解析によって地下水位を決定する。これは、デジタルツインにおいて地下水位分布がモデルの状態変数の１つであり、いくつかの地点で観測された地下水位をもとに浸透流解析を行ってモデル全体の地下水位分布を更新していくことを想定したものである。
デジタルツインでは、内部状態が更新されるのに合わせて、解析を繰り返す必要があるため、高速かつ安定な解析手法が求められる。3 次元の地震シミュレーションの高速計算手法は様々なものが開発されてきている。しかし、解析の安定性についての検討はあまりなされていない。デジタルツインのような複雑な地盤構造物モデルについて、液状化を考慮したマルチフィジックスの解析を行うと解析が不安定になってしまうことがある。そこで、破壊解析などで使用されている積分平均化理論によって空間方向に構成則の平均化を行い、さらに、高周波成分を減衰させる一般化 α 法を時間積分に用いることで、この問題を克服し、安定的に解析が実施できるような手法を開発した。地盤の構成則はある程度の体積の平均的な挙動としてモデル化されているが、数値解析上ではその構成則がポイントワイズで成立すると仮定されている。従来の解析では、分解能が低めの地盤モデルを使用していたため、この仮定をしても問題がなかったが、超高分解能な解析においてこの仮定に基づいて解析を行うと、実現象を適切に表現できなくなってしまうことがある。これによって、解析の不安定が引き起こされていると考えられる。
開発手法では構成則の平均化を行うことでこの問題を回避できるため、解析が安定化されることが期待される。
従来手法と開発手法を使って、複数の分解能で離散化した地盤モデルについて、液状化を考慮した地震応答解析を実施した。従来手法では、高分解能のモデルを使うと解析が不安定になり、収束解が得られないのに対し、開発手法では、安定に解析が実施でき、かつ解析結果が収束することが示された。開発手法を使った適用例として、大規模地震応答解析を実施した。河口付近の地形を模した 7 億自由度の地盤構造物モデルを使用し、浸透流解析によって決定された地下水位分布を使って解析を実施した。計算には AI Bridging Cloud Infrastructure (ABCI)上の NVIDIA Tesla V100 GPUs 480 台を使用した。
④-(a) B. 手法
本研究では、デジタルツイン上でのマルチフィジックス地震応答解析のための手法を開発した。デジタルツインでは、観測データなどによって内部状態が更新されていく。地震リスクを正確に評価するには、内部状態が更新されるたびに、地震応答解析を行って地震リスクを再評価する必要がある、そのため、高速かつ安定な地震解析手法が必要である。本研究では、有限要素法に基づく GPU を用いた高速な地震応答解析手法をもとに、地震応答解析の安定化手法を開発した。この手法は OpenACC を使って実装されており、ほぼすべての計算が GPU 上で実行される。また、観測データによって変化する内部状態の1 つとして地下水位分布を対象とし、地下水位分布を浸透流解析によって計算し、それを用いて液状化を考慮した地震応答解析を実施することとした。

- 89 -

1)地震応答解析の問題設定
本研究では、非構造四面体二次要素を用いた有限要素法により地震応答解析を行う。
有限要素法は複雑な geometry を表現することが可能で、非線形材料の取り扱いや、自由表面境界条件の設定も容易である。デジタルツインの地震応答解析では安定性が重要であることから、時間積分法として一般化 α 法を使用する。一般化 α 法では高周波成分に減衰がかかるため、数値誤差により発生する
高周波の波による解析の不安定化が抑制される。有限要素法と一般化 α 法により離散化された運動方程式は式(2.4.1)となる。
（2.4.1）（2.4.2）（2.4.3）

ここで、M、C、K はそれぞれ質量行列、Rayleigh 減衰行列、剛性行列であり、δu、f、q、v、a は変位
増分ベクトル、外力ベクトル、内力ベクトル、速度ベクトル、加速度ベクトルである。dt は時間ステッ
プ幅である。また、*(n)は n ステップ目の変数*を表す。αm、αf、β、γ は一般化 α 法のパラメータであり、本研究では、αm=1/3、αf=4/9、 β=(1-αm+αf)2/4、 γ=1/2-αm+αf を使用した。なお、αm=0、αf=0 とするとNewmark-β 法と一致する。
各時間ステップで、式(2.4.1)を解いて、δu を得、以下の式によって、各変数の値を更新する。
（2.4.4）（2.4.5）（2.4.6）（2.4.7）（2.4.8）
ここで、B は変位をひずみに変換する行列である。また、D は弾塑性行列、σ は全応力であり、構成則によって更新される。∫𝑉 ∗ d𝑉 は要素 e における体積積分を表す。
𝑒2)浸透流解析
デジタルツインでは、時々刻々と得られる観測データをもとにモデルの状態を更新していく。本研究では、その例の１つとして、浸透流解析により季節変動する地下水位を計算し、それぞれの地下水位で- 90 -

地震応答解析を行うことで、地震によるリスクの変動を計算することとした。従来の地震シミュレーシ
ョンでは地下水位は恣意的に、地表面から何メートルという設定をされた計算が多かった。本研究では浸透流解析によって地下水位を決定した。これにより、重要構造物のより正確な性能照査につながると期待される。浸透流解析では水を非圧縮流体と仮定し、定常流を考える。鉛直方向の流れを無視し、鉛直方向に質量保存則の式を積分することで、下記の 2 次元の支配方程式を得る。
（2.4.9）（2.4.10）（2.4.11）ただし、h は全水頭、K は透水係数テンソル、q は系外部への水の流出量である。
液状化層を透水係数が一定の不圧帯水層とし、2 次元の 1 次構造要素の有限要素法で離散化して全水頭の分布を求め、それを地下水位とした。地下水位の季節変動に対して、地震は十分短い時間で発生すると考えられるので、地震応答解析時の地下水位分布は時間変化しないものとした。

3)構成則の空間平均化による地震応答解析の安定化
デジタルツインの解析では、多数のケースについて大規模解析を行う必要があるため、解析の安定性が求められる。 しかし、デジタルツインのような複雑な地盤構造物モデルについて、液状化を考慮した
マルチフィジックスの地震応答解析を行うと解析が不安定になってしまうことがある。そこで、積分平均化理論によって構成則の平均化を行うことで、この問題を克服し、安定的に解析が実施できるような手法を開発した。
本研究では、Iai らの構成則を使用する。この構成則は、非排水条件下での地盤液状化を記述するもの
で、過剰間隙水圧モデルとマルチスプリングモデルからなっている。過剰間隙水圧モデルによって、液状化の進行度から地盤の剛性パラメータ θ(体積弾性係数、せん断剛性係数、せん断強度)が計算され、マルチスプリングモデルによって、剛性パラメータとひずみ、応力と弾塑性行列が計算される。
（2.4.12）（2.4.13）
この構成則では、塑性せん断仕事の累積値が液状化進行度のパラメータである、通常の解析手法では、
10-1m オーダーの小さな要素を使って解析を行うと、(1)液状化が進行し、(2)剛性が低下し、(3)ひずみが大きくなり、(4)塑性せん断仕事が累積され、(5)さらに液状化が進行する、というサイクルがごく一部の要素で発生し、それらの要素だけで異常に大きなひずみが生じてしまうことがある。しかし、実際に 10-1m
オーダーの局所的な領域だけで液状化進行することは考えづらい。問題は、従来の解析では構成則がポイントワイズで成立すると仮定していることにある。そもそも、地盤は不均質であり、構成則はある一定の体積における平均的な挙動をモデル化したものである。通常の解析手法では 1m～10m オーダーの要
素サイズで解析をしていたため、構成則がポイントワイズで成立すると仮定したとしても、そこまで大- 91 -

きな問題にはならない。しかし、デジタルツインのような高分解能な解析では、構成則の適用範囲外の短いレングスケールで構成則が計算されてしまい、解析の不安定化につながっていると考えられる。本
研究では、構成則をポイントワイズに成立するものではなく、ある体積における平均的な応力ひずみ関係であるととらえることで、この問題の克服する手法を開発した。構成則の平均化は以下の手順によって行う。
1. 従来の構成則と同様に、過剰間隙水圧モデルによって、各位置 x における剛性パラメータ θ(x)を計算する。
̅ (𝒙)を以下のように計算する。
2. 積分平均化理論によって、剛性のパラメータの非局所化量𝜽
（2.4.14）
ただし、重み関数 α(x、ξ)は x=ξ の時に最大値をとり、r=||x-ξ||について単調減少な関数で、を満たす。本研究では、以下の関数を使用する。
（2.4.15）（2.4.16）

ここで、R は平均化の範囲を表すパラメータで、本研究では、R=1m を使用した。
3. 剛性パラメータの代わりに、剛性パラメータの非局所化量を使って、応力と弾塑性行列を計算する。
剛性パラメータ以外は、オリジナルの構成則のマルチスプリングモデルによる計算と同じである。
（2.4.17）（2.4.18）
この手法により、局所的に一部の要素で液状化により塑性せん断仕事が累積されたとしても、剛性の低下が抑えられる。そのため、前述のひずみ集中が発生するサイクル((1)液状化が進行し、(2)剛性が低下し、(3)ひずみが大きくなり、(4)塑性せん断仕事が累積され、(5)さらに液状化が進行する)において、(1)から(2)が起こりにくくなり、解析の不安定化が抑制されると期待される。
④-(a) C. 数値実験
本研究で開発した安定化手法の有効性を示すために、数値実験を行った。まず、従来の手法では不安定化してしまう高分解能な地盤モデルを使った解析が、本研究の安定化手法を用いることで安定的に解析でき、かつ解析結果が収束することを示す。次に、適用例として、河口付近の地形を模した地盤構造物モデルを用いた地震応答解析を行った。浸透流解析を行って地下水位分布を計算し、その分布をもとに地震応答解析を実施した。全ての解析において、1995 年兵庫県南部地震の際に観測された地震波を入力波として使用した。

- 92 -

1)従来法との収束性の比較開発手法により、高分解能なモデルでの地震応答解析が安定的に実施できるようになったことを示す。
図 2.4.1 に示す地盤モデルを使用し、液状化を考慮した地震応答解析を行った。
物性パラメータは表 2.4.1 に示す。地下水位は地表面から 2m に設定した。時間ステップ幅 dt=0.001sである。空間離散化の分解能が 1m、0.5m、0.25m、0.125m の 4 ケースについて、従来手法と開発手法で解析を行い、解の収束性を確認した。ここで、従来手法と開発手法の違いは構成則による応力と弾塑性行列の計算方法と時間積分法である。従来手法では、空間平均化していない剛性パラメータを使って、式(2.4.12、2.4.13)により、応力と弾塑性行列を計算し、時間積分法に Newmark-β 法(式(2.4.3)で、αm=0、αf=0、
β=1/4、 γ=1/2 としたもの)を使用した。開発手法では剛性パラメータの非局所量を使って、式(2.4.17、
2.4.18)により、応力と弾塑性行列を計算し、時間積分法に一般化 α 法を使用した。

図 2.4.1 収束性確認に使用された地盤モデル
表 2.4.1 Model Properties

図 2.4.2 に地表面中心部における変位の時刻歴を示す。従来法では、分解能を小さくしても、解析結果が収束しない。また、分解能 0.25m、0.125m のケースでは解析が t=5.6s 付近で異常終了した。図 2.4.3 に- 93 -

分解能 0.125m での解析における解析での異常終了時の変位分布を示す。少数の要素だけが異常な変形をしており、ひずみが集中してしまっている。一方、開発手法を使った場合、分解能 0.125m のケースでもひずみが集中することなく解析が実行できており、解析結果も収束している。開発手法により、高分解能なモデルを使っても安定に解析が実行でき、かつ、収束解が得られることが示された。

図 2.4.2 地表面中央部での変位の時刻歴

図 2.4.3 解像度 0.125ｍでの解析における異常終了時の変位ノルムと変形

2)適用例：高分解能地盤構造物モデルによる地震応答解析
河口付近の地形を模した高分解能な地盤構造物モデルを用いた地震応答解析を行った。分解能は 0.5mであり、構成則の平均化を行わない従来法を用いて計算した場合、解析が不安定になり得る分解能であ- 94 -

る。使用した地盤構造物モデルを図 2.4.4 に示す。モデルの自由度は 702,192,969、要素数は 172,929,616である。モデルの物性を表 2.4.2 に示す。液状化特性のパラメータは表 2.4.1(c)に示されたものを使用した。
液状化層を不圧帯水層とみなして地下水の浸透流解析を行い、その解析結果を地下水位分布に設定した。
地下水位の分布を図 2.4.5 に示す。開発手法を使用し、時間ステップ幅 dt=0.001 s で 20 秒間分(20,000 time
steps)の液状化を考慮した地震応答解析を行った。480 MPI processes x 1 GPU per 1 MPI process = 480 GPUs
による並列計算によって計算した。計算は、GPU を搭載したスーパーコンピュータ AI Bridging Cloud
Infrastructure (ABCI)の 120 計算ノードを使用した。ABCI の各計算ノードは NVIDIA Tesla V100 GPU 4 台と Intel Xeon Gold 6148 CPU 2 台を搭載している。計算時間は 4 時間 20 分だった。地表面における残留変位の分布を図 2.4.6 に示す。開発手法によって、大規模かつ複雑な地盤構造物モデルにおける高分解能な地震応答解析を安定的に実施できるようになったことが示されている。

図 2.4.4 Model used in the application
表 2.4.2 Model properties of the application- 95 -

図 2.4.5 浸透流解析により計算された地下水位分布

図 2.4.6 適用例おける t =20s での地表面における変位④-(a) D. まとめ
本研究では、デジタルツインのための地震応答解析手法を開発した。デジタルツインでの解析は高速であるだけでなく、安定であることも求められる。しかし、液状化などのマルチフィジックスを考慮した解析では、デジタルツインで使用される超高分解能な地盤構造物モデルを使うと、不安定化してしま
うことがしばしばある。本研究では、開発された大規模数値解析コード群の実展開を踏まえた解析手法の安定化及び有効性の確認を目指して、構成則の平均化による安定化手法を開発した。この手法によって、超高分解能なモデルの解析においても、安定的に解析が実施でき、かつ収束解を得ることができるようになった。適用例として、7 億自由度の地盤構造物モデルについての地震応答解析を実施した。より正確な地震リスク評価を目指して、地下水位を浸透流解析によって計算し、その地下水位分布に基づいて地震応答解析を実施した。本研究では、特定の液状化構成則についての実装を行ったが、開発手法は構成則とは独立のものであり、他の地盤構成則にも適用できる汎用的なものである、本手法による解析
結果の信頼性向上のための実観測データとの融合は今後の課題である。また、3 次元での非定常浸透流解析や、排水条件下での液状化解析などを実施することにより、さらなる信頼性の向上が期待できる。
④-(b) GPU 上での液状化解析の高速化と解析結果を用いたサロゲートモデルの構築④-(b) A. はじめに
Physics-based の地震シミュレーションは地震被害の推定や軽減に資すると期待されているが、特に 3- 96 -

次元の解析では計算コストが膨大である。このコストを克服するために、我々は、有限要素法を用いた大規模な 3 次元地震応答解析のための様々な手法を開発してきた。本研究では、マルチフィジックスである液状化現象を考慮した波動伝播解析を対象として、開発された大規模数値解析コード群の実展開を踏まえた解析手法の安定化及び有効性の確認を行った。具体的には、マルチフィジックス問題の特徴を考慮した計算負荷分散手法によって高速化を図る（開発手法は CPU のみを使用した手法と比較して 10倍以上の高速化を達成し、解析に必要な計算資源を大幅に削減した）とともに、安定的に解析が可能となるようにした。さらに、安定的かつ高速に解析を可能と出来るように、解析結果を用いた機械学習に
より地震シミュレーションのサロゲートモデルを構築することにより、より短時間で安定的に地震応答の評価が可能になることを示した。
④-(b) B. GPU を使用した高速な地震シミュレーション
本研究の支配方程式は、以下の地盤の運動方程式(2.4.19)である。
（2.4.19）
ただし、σ は応力、u は土骨格の変位、ρ は地盤の密度、f は外力、t は時刻である。また、地盤の構成則として、Iai のマルチスプリングモデルを使用した。式(2.4.19)は有限要素法と Newmark-β 法によって離散化され、100 万～10 億自由度の疎行列方程式となる。時間ステップごとに係数行列が更新されるこの方程式を数万時間ステップにわたって解く必要がある。開発手法では可変的前処理付き共役勾配によって式(2.4.19)を解いている。前処理において、マルチグリッド法と精度混合演算を用いることで計算コストを削減している。解析全体は、MPI と OpenACC により並列化されている。マルチフィジックス解析の特徴を考慮した領域分割の改良と、要素リオーダリングを採用し、並列効率が高く、かつ GPU アー
キテクチャに適したアルゴリズムを開発した。また、21 ビットの低精度浮動小数点数 FP21 を前処理計
算における MPI 通信に使うことでさらなる高速化を達成した。対象問題はマルチフィジックス現象を考慮した複雑な問題あり、解析は不安定になり得るが、Section④-(a)の手法により安定した解析が可能になっている。
スーパーコンピュータ AI Bridging Cloud Infrastructure (ABCI) の 1 計算ノード(NVIDIA V100 GPUs を 4台搭載) 上での性能計測において、開発手法は CPU のみを使用した手法と比較して 10 倍以上の高速化を達成した。開発手法により、89,146,716 自由度の地盤構造物モデルについての 30,000 時間ステップの解析を 13 計算ノード(ABCI)を使って 3 時間 33 分で実施できるようになった。
この解析を Intel Xeon Phi 7250
(Knights Landing) CPU を搭載したスーパーコンピュータで GPU を使用せずに実施した場合、128 計算ノードを使って 14 時間 37 分かかった。
④-(b) C. 解析結果を学習に使用したサロゲートモデルの構築
解析が安定化・高速化したことにより、より多くの解析データが得られるようになった。このデータを用いて、機械学習を行い、地震シミュレーションのサロゲートモデルを構築した。これにより、さらに高速な安定した地震応答の評価を可能にした。サロゲートモデルは 2-4 層の全結合層からなるニューラ
ルネットワークである。液状化層の厚さや地表面からの深さが液状化のしやすさに影響するという経験- 97 -

則に基づき、各層の厚さを入力データとし、その地点で液状化する(液状化の指標の１つである過剰間隙
水圧比が 0.9 を超える)かどうかを出力するニューラルネットワークを構築した。ラテン超方格法により対象とする大規模地盤モデル(0.9 億自由度)の層厚データのサンプリングを行い、その層厚データを使用
して 10 個の小規模地盤モデル(100 万自由度)を作成した。この小規模モデルの解析結果を学習データに使用し、機械学習を行った。これにより、対象の地盤モデルの統計量をもとに学習データが生成され、少ないサンプル数で効率的な学習が実現できた。図 2.4.7 に推定結果を示す。９割以上の正答率を実現しており、サロゲートモデルの有効性が示されている。

図 2.4.7 物理シミュレーションとサロゲートモデルの解析結果の比較④-(b) D. まとめ
本研究では、GPU を使って大規模な液状化解析の高速化・安定化を実現した。本研究では、1 つの構
成則についての実装が示されているが、他の構成則にも適用可能な汎用的な手法である。また、多数回の解析を高速・安定に実施できるようになり、大量の解析データを得ることが出来るようになったことを活かして、この解析データをもとに作成したサロゲートモデル構築により、さらに高速に安定的に地震応答の評価が可能になった。これにより、大規模領域について、通常の Physics-based の地震シミュレーションでは実現困難な数のシナリオに対する地震応答の評価が可能になると期待される。
（３）開発した計算アプリケーション群を実務で活用するための環境整備（１）において整備している国の被害想定のためのハザード評価に使われるものと同等の長周期地震
動計算アプリケーション群の実務への活用を進めるため、連携機関の各企業が独自の問題設定で実施するためのサポート体制を、海洋研究開発機構が中心となって整備している。対象としているのは計算ア
プリケーション群の基幹数値解析コンポーネントである E-wave FEM の計算環境である。本年度は、昨年
度に構築した計算環境を改良し、より現実的な地下構造を対象とした問題への適用をしやすくした。更に、当該計算環境を有効に利用して頂くため、連携企業メンバーを対象に、ウェビナー形式（Web 上のオンラインセミナー）のチュートリアルを実施した。その上で、連携企業メンバーから E-wave FEM 計算
環境に関するフィードバックを引き続き受け、計算環境の改良をさらに進めた。本稿①にて、これらの取り組みについて述べる。また、連携企業では、企業でニーズのある地震動計算への活用を進められている。その中で得られた成果について②で述べる。
①実施した計算環境整備と連携企業サポートの概要- 98 -

背景として昨年度の取り組みの概要を述べた後に、本年度の取り組みについて説明する。
①-(a) 令和 2 年度までの取り組みの概要
令和 2 年度までに、計算アプリケーション群の基幹数値解析コンポーネントである E-wave FEM の計算環境として、大規模共有メモリ計算機からなる計算システム（以下、LM システムと記載 ）及び高並列
スーパーコンピュータからなる計算システム（以下、SC システムと記載）を連接する計算システムを構築した。LM システムを構成する大規模共有メモリ計算機として、海洋研究開発機構が所有し理化学研究
所計算科学研究センターに設置された大規模共有メモリ計算機”Jasper”を使用した。また、SC システムを構成するスーパーコンピュータとして、東京大学情報基盤センターと筑波大学計算科学研究センターが共同運営する、
最先端共同 HPC 基盤施設
（JCAHPC：Joint Center for Advanced High Performance Computing）
の

共

同

利

用

（

高

並

列

）

ス

ー

パ

ー

コ

ン

ピ

ュ

ー

タ

”Oakforest-PACS”(https://www.cc.u-tokyo.ac.jp/supercomputer/ofp/system.php)を使用した。実際の E-wave

FEM の計算においては、LM システムにプリプロセスを、SC システムに E-wave FEM による地震動計算（以下、本解析）、及びポストプロセスを割り当てている。尚、上記プリプロセス、本解析、及びポストプロセス各々は、更に細かな段階的解析処理ステップ（以下、処理ステップ）に分割されている。
連携企業のユーザに対し、E-wave FEM に関する 5 回のウェビナー形式のチュートリアルを提供した。
E-wave FEM によって、比較的単純な条件下（簡単化した地盤モデルに対して点震源あるいは地盤モデルの境界面に設定された面震源から発生する）での長周期地震動計算及び解析について説明した。
チュートリアルでは、計算システム全体を説明すると共に、ハンズオン（遠隔実習）により、実際に大規模共有メモリ計算機”Jasper”、及び高並列スーパーコンピュータ”Oakforest-PACS”に SSH 接続し、プリプロセス、本解析、ポストプロセスの一連の流れを経験して頂いた。またユーザからは、一連のチュートリアル及び E-wave FEM 利用に対するフィードバックを受けた。
①-(b) 令和 3 年度の取り組み
令和 3 年度は、令和 2 年度の成果を踏まえて、計算環境の拡充に伴うシステム変更、及びユーザからのフィードバックを踏まえたアプリケーション群の改良、機能拡充を実施した。また、連携企業サポー
トとして、連携企業ユーザに対して 6 回目のチュートリアルを実施した（令和 3 年 10 月実施）。このチュートリアルでは、計算環境の拡充によるシステム変更、及びアプリケーション群の機能拡充に伴う利用手続きの変更を把握して頂くことを主な目的とした。チュートリアル後から現在まで、ユーザからの更なるフィードバックを受け、アプリケーション群に対する機能追加、機能改良、バグ修正を継続的に行なっている。以下に、本年度実施した 1) 計算環境の拡充に対応したシステム変更、2) アプリケーショ
ン群の改良による機能拡充、3) チュートリアルの概要、4) チュートリアル後の更なるアプリケーション
群の改良の 4 点について順に述べる。
1) 計算環境の拡充によるシステム変更
新たな計算環境として、「富岳」（理化学研究所）の本格運用開始（R3 年 3 月）、地球シミュレータ第 4
世代 (ES4:海洋研究開発機構: https://www.jamstec.go.jp/es/jp/system)の GPU(Graphical Processing Units)ノードの運用開始（R3 年 6 月）、と利用可能計算機が拡充する環境となった。これら計算機利用条件の変更を受け、従来の計算環境に複数の計算機を追加し、利用可能とするシステム変更を行なった。具体的には、LM システムを構成する大規模共有メモリ計算機として、ES4 の GPU ノードを利用可能とした。また、- 99 -

SC システムの高並列スーパーコンピュータとして、「富岳」及び上記 ES4 の VE（Vector Engine）ノードを利用可能とした。その結果として、プリプロセス（LM システム利用）に Jasper 及び ES4（GPU）、本
解析とポストプロセス（SC システム利用）に Oakforest-PACKS、「富岳」及び ES4（VE）を利用できるよ
うになった。また、上記システムの変更に伴い、各処理ステップのコントロール（パラメータ設定）ファイル及びその格納ディレクトリの整理を行なった。具体的には、変更対象となるコントロールファイ
ルを利用する計算機毎に作成し、その適用プロセスに応じて LM システム専用、SC システム専用、LM/SC
システム共用のディレクトリに各々格納した。この整理によって、ユーザがこれまでと異なる計算機の組み合わせを利用する場合であっても、見通し良くパラメータの設定をすることが期待できる。
2) アプリケーション群の機能拡充
アプリケーション群の機能拡充に関して、拡充する 3 つの機能、それらに伴う新規パラメータの追加、及び手続きの変更の順に説明する。まず、拡充する 3 つの機能とは、第 1 に実地盤データの導入機能、第 2 にメッシュモデル（地盤データの情報を有限要素節点に割り当てた有限要素モデル）外周部（底面を含む）内側への吸収帯の設定機能、第 3 に強化されたチェック機能である。
第 1 の実地盤データの導入機能とは、新たに設定した地盤層データファイルを介して、外部から取り入れた実地盤データ（層境界形状と層内の物性値）を本解析用の有限要素節点に割り当て、本解析の際に実地盤データを反映させるものである。実地盤データは、例えば、地震動計算に広く利用されている
数値標高データ（DEM データ:https://fgd.gsi.go.jp/download/ref_dem.html）を利用できることを想定している。そのため、地理座標系の DEM データファイルと E-wave FEM システムが直接読み込む地盤層データ
ファイル（FEM 計算用直交座標系）との間で、ファイル変換をするツールも併せて開発した。DEM デー
タに限らず何れの実地盤データであっても、実地盤データ情報を地盤層データファイルに移すことができれば本解析に実地盤データを導入できる。
第 2 のメッシュモデル外周部内側に吸収帯[1]を付加できる機能は、本解析の際にメッシュモデル外周部からの反射波の影響を低減させる目的で導入した。吸収帯厚及びその物性は、予め定められた吸収帯に関わる各種パラメータから自動的に計算され、ユーザの負担なしに吸収帯を設定できる仕組みとした。
吸収帯の付加の有無はユーザが簡単にパラメータ（on/off のフラッグ）で選択できる。また、吸収帯厚の調整を手動で行なう機能も追加した（後述）。
第 3 の強化されたチェック機能に関しては、機能拡充に伴うモデル変更部分をチェックする 2 つの機能を追加した。一つ目の機能は、メッシュモデル作成時の有限要素形状の品質チェック機能である。実
地盤データを導入する場合、その層境界形状によっては、層境界近傍の有限要素形状の品質が計算に不適切（要素の体積が 0、アスペクト比が過大、要素同士の不整合が発生する等）になる場合がある。その品質を示すインデックスをファイルに出力し、不適切性を十分に改善できない場合はメッシュモデルの設定をやりなおすことをユーザに促すシステムとした。二つ目のチェック機能として、メッシュモデルと吸収帯との位置関係を告知する機能を設けた。吸収帯設定はメッシュモデル作成と同時に行われ、その際にメッシュモデルサイズと吸収帯との位置関係を標準出力に示すことで、メッシュモデルに対し吸収帯配置が適切であるかをユーザが把握しやすい形とした。
（例えば、震源域或いは重要な観測点に相当する計算節点が吸収帯に含まれるか否かを判断できる。）この表示により、メッシュモデル作成の初期段
階で計算領域（メッシュモデル全体のサイズ）を取り直すこともできる。また吸収帯厚は、他の条件が揃った本解析直前であっても（メッシュモデルを変更せずに）変更できる設定も追加した。

- 100 -

上記 3 つの機能拡充（実地盤データ導入、吸収帯付加、チェック機能強化）に伴い、関連する新たなパラメータを設定し、対応する処理ステップのコントロールファイルに追加した。計算手順の変更として、機能拡充に伴い、プリプロセス、及び本解析＋ポストプロセス内の各処理ステップを再構成し、各処理ステップの実行手順の変更が発生した。特に、実地盤データの読み込みから、計算用のメッシュモ
デル作成までの手順は、実行内容毎（DEM データから読み込み可能な地盤層データファイルへの変換、
地盤層データファイルからのメッシュモデル作成、メッシュモデルの品質チェック、品質内容に応じた
メッシュモデルの局所再分割、再度のメッシュモデルの品質チェック）に細かく処理ステップに分割した。その結果として、ユーザが各処理ステップにおいて、処理内容の理解、確認、及び誤りの発見が容易になるように配慮した。また、メッシュモデルとその物性値の 3 次元表示、及び FEM 計算結果の 2 次
元・3 次元波形表示は、前者はメッシュモデル作成の終わったプリプロセス内の最終処理ステップで、後者はポストプロセス内の最終処理ステップで実施する手順へと整理、変更した。2 つの可視化（プリプロ
セス最終処理ステップでのメッシュモデル可視化、ポストプロセス最終処理ステップでの計算結果の可視化）の分離は、ユーザの数値実験デザインの再検討、或いは実験設定の誤り発見等を容易にすることを意図している。
以上の 1) 計算環境の拡充に対応したシステム変更、2) アプリケーション群の改良による機能拡充に関しては、実装が完了した段階で、その変更内容を基本マニュアルに反映している。
3) 第 6 回チュートリアルの実施
1) 2)の計算環境更新に伴い、令和 3 年 10 月 5 日に、本課題の連携企業ユーザに対して 6 回目のチュー
トリアルを開催した。チュートリアルの内容の概要を表 3.1.1 に示す。チュートリアルでは、前回チュー
トリアルからおよそ一年が経過していることを踏まえ、E-wave FEM 使用方法の根幹部分（有限要素計算
用メッシュの生成、並列計算用の領域分割、スーパーコンピュータ上での高並列地震動計算、計算結果の可視化まで）をウェビナー形式で説明した。更に、それまでに実装した計算環境の拡充に伴うシステ
ム変更、及びアプリケーション群の機能拡充に伴う利用手続きの変更部分に焦点を当てて、説明を加えた。

- 101 -

表 3.1.1 チュートリアル実施内容
日付、形式
第 6 回チュートリアル
令和 3 年 10 月 5 日、ウェビナー形式

内容（E-wave FEM 計算手続きと新規事項の説明）
1. 新規計算環境の紹介、各計算機との接続説明
2. チュートリアル用問題の紹介
3. プリプロセス
・地盤モデル設定
・地盤物性設定
・地盤層データファイル等の作成（新規）・メッシュモデルの作成、吸収帯の付加（新規）・メッシュモデルの品質チェック（新規）・メッシュモデルの細分化
・メッシュモデルの分割（並列計算用）
・分割済みメッシュモデルをスーパーコンピュータに転送
・メッシュモデルの可視化
4. 本解析
・解析、震源及び観測点の設定
・観測点等の最近傍有限要素節点の抽出
・本解析の実⾏
5. ポストプロセス
・ポストプロセスの実⾏・計算波形の 2 次元・3 次元表示
・計算結果の各種フォーマットへの変換、ファイル出力4) チュートリアル後の更なるアプリケーション群の改良
第 6 回チュートリアルに関するユーザからのフィードバックを受けて、更に機能の追加、改良、バグ修正を実施した。機能の追加としては、以下 4 種類の機能があげられる。第 1 に、地盤層データファイ
ル作成するステップにおいて、入力する DEM データのチェック機能を追加した。第 2 に、メッシュモデ
ル生成するステップにおいて、ソフトウェア上でのモデルサイズ制限を緩和する機能を追加した。第 3に、メッシュを細分化するステップにおいて、SC システムでの FEM 計算使用ノード数の目安を提案する機能を追加した。メッシュ細分化ステップで確定する有限要素モデルの要素数に応じて、FEM 計算におけるメモリ不足を回避できると思われる最小ノード数の経験的目安を、使用ノード数下限として提案する。第 4 に、多数の点震源で震源域を構成する場合、点震源の最近傍有限要素節点を探索する処理を高速化するためのスレッド並列化を可能とする機能を追加した。機能の改良としては、DEM データから
メッシュ生成を行うステップでの、地盤構造データ補間アルゴリズムの改良がある。これは、特定の条件下で生成されたメッシュに不自然な形状が生じる場合があることがわかっていたため、この現象を抑制するために導入したものである。バグ修正に関しては、地表面波動場可視化用データ（vtk ファイル）
作成時のバグ修正、面震源作成プログラムのバグ修正、モデル境界での吸収帯幅設定のバグ修正、を実施している。
これらの機能追加、改良、バグ修正によるアプリケーション群の高機能化、安定化には、連携企業ユーザからのフィードバックが極めて重要であり、今後とも連携企メンバーとの連携を重視していく。

- 102 -

② 企業においてニーズのある地震動計算への適用―複数の発生シナリオを考慮した長周期地震動評価―
チュートリアルにより E-wave FEM の使い方を習得された連携企業メンバーによって、企業でニーズのある地震動計算への活用が進められている。その中で本年度得られた成果について述べる。
今後発生が懸念される千島海溝沿いの巨大地震では、内閣府による被害想定結果[2]が示され、甚大な津波被害が予測されている。一方で長周期地震動の発生により、札幌の超高層建築物や免震建築物、苫小牧の石油タンク等の固有周期が長い建築物への影響も大きいと考えられ、当該地域での長周期地震動の評価が重要になっている。また、震源特性は必ずしも一意ではなく不確実性を有しており、その選択により予測結果に大きなばらつきが生じるため、多様なシナリオを考慮した検討が必要である。そこで、千島海溝沿いの巨大地震を対象に、「富岳」を用いて 3 次元有限要素法により長周期地震動を評価した。
また、震源特性の不確実性を考慮した複数の震源モデルを構築し、ばらつきの大きさや各震源特性が長周期地震動に与える影響を分析した。
地震動評価には、地震動計算用大規模有限要素計算プログラム「E-wave FEM」[3][4][5]を用いた。計算
資源は、地震動評価には「富岳」および「Oakforest-PACS」を、地盤メッシュ分割には「地球シミュレー
タ（ES4）」を用いた。
地盤モデルは、千島海溝沿いの長さ 1,000km、幅 420km、深さ 150km の範囲とし、地下構造モデルには地震本部の全国 1 次地下構造モデル[6]を用いた。なお、全国 1 次地下構造モデルの範囲外となる東経147°以東は海溝軸に沿って外挿した。対象周期 3.3 秒以上、1 波長 5 要素以上で要素分割を行い、約 14.7
億節点、10.8 億要素となった。図 3.2.1 に地盤モデルのプレート境界深度、Vs1,500m/s 層下面深度を、表
3.2.1 に地盤物性を示す。震源モデルは、十勝沖から択捉島沖を震源域とする地震本部[7]の「超巨大地震（17 世紀型）」を想定し、表 3.2.2 に示す震源特性に不確実性を考慮した 28 ケースを設定した。図 3.2.1に設定した震源モデルを示す。震源パラメータは、内閣府の被害想定[2]を参考にレシピ[8]に従って設定した。表 3.2.3 に震源パラメータを示す。また、山本・吉村[9]と同様に、すべり速度時間関数は 2 つの三角形で表現した。図 3.2.2 に、すべり速度時間関数の模式図を示す。ここで、地盤モデルに外挿範囲を含むことから、妥当性の検証として 2001 年 5 月 25 日択捉島南東沖の地震について、観測記録と解析結果を比較した。F-net のメカニズム解、気象庁の震源位置を用い、点震源でモデル化した。図 3.2.1 に震源位置を、表 3.2.4 に各地震の諸元を、図 3.2.3 に比較結果を示す。堆積層の厚さが異なる 2 地点とも振幅を概ね再現することから、地震動評価に使用できると判断した。
図 3.2.4 に、根室、苫小牧、札幌について擬似速度応答スペクトル（pSv）を示す。合わせて 2003 年十勝沖地震の観測記録を示す。一部ケースでは 2003 年十勝沖地震や告示 Lv.2（構造物の耐震設計に用いる
入力地震動で，現在から将来にわたって当該地点で考えられる最大級の強さをもつ地震動）を上回る地
震動レベルとなった。図 3.2.5 に、pSv が+1σ に最も近いケースの加速度波形を示す。継続時間が長く、免制震装置や鉄骨の繰り返し変形による影響が大きくなる可能性がある。図 3.2.6 に、pSv のばらつき（自然対数標準偏差）の大きさを示す。0.4～1.2 の範囲となり、熊本地震前震[10]や相模トラフ[11]を対象とした既往検討と比較して、ばらつきが大きい傾向がみられる。本検討は震源域が広いため、震源特性のうち破壊開始点、アスペリティ位置といった幾何学的な特性の影響がより大きくなっている可能性があ- 103 -

る。図 3.2.7、図 3.2.8 に、地盤モデルを外挿していない範囲について、平均 pSv、ばらつきの大きさの空間分布を示す。pSv は堆積層が厚い地域で大きく、深部地盤の影響がみられる。震源域が広い場合、遠方の震源からの地震波は斜め入射となり、表面波が励起されやすく、地盤の影響が大きくなる可能性がある。一方で、ばらつきの大きさは pSv の大小と必ずしも対応していないことから、震源特性の違いが影響している可能性が考えられる。
ここで、不確実性を考慮した各震源特性が地震動強さに与える影響を確認するため、一つの震源特性を固定した各ケースの pSv の平均値と、全ケースの pSv の平均値の比率を算出し、各震源特性について比較した。図 3.2.9 に比率を示す。比率の常用対数の絶対値（図中 1.0 からの距離）が地震動強さへの影
響の大きさであり、正負は影響の正負に対応する。破壊開始点、アスペリティ位置は、地盤の影響を含み大きく評価される可能性があるため、今後深部地盤も含めた分析が必要であるが、本検討で不確実性を考慮した特性では、各地点とも破壊開始点の影響が大きい。また、アスペリティが浅い位置にある方が、地震動強さが大きくなる傾向がみられた。
本検討では、千島海溝沿いの巨大地震を対象として、「富岳」を用いて地震動計算用大規模有限要素計
算プログラム「E-wave FEM」により長周期地震動を評価した。その際、内閣府の被害想定、地震本部の長期評価を参考に、震源特性の不確実性を考慮した 28 ケースの震源モデルを構築し、ばらつきの大きさ
や各震源特性が長周期地震動に与える影響を分析した。長周期地震動評価の結果、一部のケースでは告示 Lv.2 や、2003 年十勝沖地震の観測記録を上回る地震動となる可能性があることを確認した。また、計算された加速度時刻歴波形から、長周期地震動の継続時間が長く、免制震装置や鉄骨の繰り返し変形による影響が懸念される。長周期地震動の強さは、堆積層の厚さが影響するが、震源特性では破壊開始点の影響が大きい。
「富岳」を用いた大規模な長周期地震動計算が行えるようになったことにより、このよ
うな長周期地震動評価の信頼性を向上させるこのような取り組みが可能となったといえる。

図 3.2.1 地盤モデル・震源モデル- 104 -

表 3.2.1 23 層の地盤物性
Layer

Vp[km/s]Vs[km/s]ρ[g/cm3]

No.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

Qualityfactor

1.7
1.8
2.0
2.1
2.2
2.3
2.4
2.7
3.0
3.2
3.5
4.2
5.0
5.5
5.8
6.4
7.5
5.0
6.8
8.0
5.4
6.5
8.1

0.35
0.5
0.6
0.7
0.8
0.9
1.0
1.3
1.5
1.7
2.0
2.4
2.9
3.2
3.4
3.8
4.5
2.9
4.0
4.7
2.8
3.5
4.6

1.80
1.95
2.00
2.05
2.07
2.10
2.15
2.20
2.25
2.30
2.35
2.45
2.60
2.65
2.70
2.80
3.20
2.40
2.90
3.20
2.60
2.80
3.40

70
100
120
140
160
180
200
260
300
340
400
400
400
400
400
400
500
200
300
500
200
300
500

表 3.2.2 考慮した震源特性の不確実性
検討ケース

震源特性

ケース数

アスペリティサイズ

2 ケース

大きい（実線）・小さい（点線）

アスペリティ位置

2 ケース

浅い（緑線）・深い（青線）

破壊開始点位置7 ケース

震源域端部・各震源域中央（星印）（カッコ内線種・記号は図 3.2.1 に対応）

- 105 -

表 3.2.3 震源パラメータ
アスペリティ小さいケース
震源域

S

M0
2

[km ]

D

[Nm]

[m]

τ
[s]

アスペリティ大きいケース
S
2

[km ]

M0

D

τ

[Nm]

[m]

[s]

1

2,413 1.92E+21 19.4 8.98

4,969 3.95E+21 19.4 12.88

色丹島沖

2

2,050 1.50E+21 17.8 8.27

4,188 3.06E+21 17.8 11.83

及び3

763 3.40E+20 10.9 5.05

1,569 7.01E+20 10.9 7.24

択捉島沖

4

775 3.49E+20 11.0 5.09

1,575 7.05E+20 10.9 7.25

背景
5
根室沖

十勝沖

全域

67,819 2.11E+22

7.6 26.24

2,131 9.16E+20 10.5 8.44

61,519 1.68E+22

6.7 26.24

4,369 1.88E+21 10.5 12.08

6775 2.01E+20

6.3 5.09

1,569 4.04E+20

6.3 7.24

背景20,413 3.37E+21

4.0 31.57

17,381 2.20E+21

3.1 31.57

7

700 2.43E+20

8.5 4.84

1,469 5.13E+20

8.5 7.00

8

1,919 1.10E+21 14.0 8.01

3,975 2.28E+21 14.0 11.52

9781 2.87E+20

8.9 5.11

1,613 5.90E+20

8.9 7.34

背景32,969 7.10E+21

5.3 37.87

29,313 5.34E+21

4.4 37.87

アスペリティ

12,306 6.85E+21

25,294 1.41E+22

背景121,200 3.16E+22

108,213 2.44E+22

合計

133,506 3.85E+22

133,506 3.85E+22

図 3.2.2 すべり速度時間関数の例（山本・吉村[9]より引用）- 106 -

表 3.2.4 択捉島南東沖の地震諸元
諸元

値

日時

2001/05/25 09:40:54

緯度44.3185

経度

148.7072

深さ [km]30

走向 [deg]227

傾斜角 [deg]25

すべり角 [deg]108

M0 [Nm]

8.71E+18

Mw

6.6

図 3.2.3 観測と解析の比較結果

図 3.2.4 擬似速度応答スペクトル（EW）- 107 -

図 3.2.5 加速度時刻歴波形の結果（EW）

図 3.2.6 ばらつきの大きさ

図 3.2.7 擬似速度応答スペクトルの空間分布（5 秒、EW）- 108 -

図 3.2.8 ばらつきの大きさの空間分布（5 秒、EW）

図 3.2.9 震源特性の地震動強さに与える影響（EW）[1] Cerjan, C., Kosloff, D., Kosloff, R., & Reshef, M. (1985). Discrete acoustic and elastic wave equations.
Geophysics, 50(4), 705–708. https://doi.org/10.1190/1.1441945[2] 内閣府：日本海溝・千島海溝沿いの巨大地震の被害想定について, 2021.
[3] T. Ichimura, M. Hori, J. Bielak A hybrid multiresolution meshing technique for finite element three
dimensional earthquake ground motion modelling in basins including topography, Geophysical JournalInternational, Volume 177, Issue 3, Pages 1221-1232, June 2009[4] T. Ichimura, K. Fujita, S. Tanaka, M. Hori, M. Lalith, Y. Shizawa, and H. Kobayashi.：Physics-based urban
earthquake simulation enhanced by 10.7 BlnDOF x 30 K time-step unstructured FE non-linear seismic wave
simulation. In SC’14 : Proceedings of the International Conference on High Performance Computing,
Networking, Storage and Analysis, pp 15–26、2014.
[5] Ichimura, T., Fujita, K., Quinay, P. E. B., Maddegedara, L., Hori, M., Tanaka, S., Shizawa Y., Kobayashi H.
and Minami, K. : Implicit nonlinear wave simulation with 1.08T DOF and 0.270T unstructured finite elements- 109 -

to enhance comprehensive earthquake simulation. In SC'15: Proceedings of the International Conference forHigh Performance Computing, Networking, Storage and Analysis, pp.1-12, IEEE, 2015.
[6] 地震本部：全国１次地下構造モデル(暫定版), 2012.
[7] 地震本部：千島海溝沿いの地震活動の長期評価（第三版）, 2017.
[8] 地震調査研究推進本部：震源断層を特定した地震の強震動予測手法 (「レシピ」)、2020[9] 山本優・吉村智昭：3 次元大規模 FEM による東海・東南海・南海連動地震の長周期地震動シミュレーション, 日本建築学会構造系論文集, No.77, Vol.677, pp.1055-1064, 2012.
[10] 西本昌, 山本優, 縣亮一郎, 藤田航平, 市村強, 堀高峰：長周期地震動に対する不確定性を考慮した断層
モデル地震波作成方法の妥当性検討、日本建築学会大会(東海)学術講演梗概集, 構造 2, pp. 397-398,
2021.
[11] 山本優, 西本昌, 五十嵐さやか, 内山泰生, 糸井達哉：相模トラフ沿いの地震のばらつきを考慮した地震動解析とリスク評価（その 1）計算概要と検証、日本建築学会大会(北海道)学術講演梗概集,
2022.(投稿中)（４）プロジェクトの総合的推進
プロジェクトの総合的推進では、1) 参画各機関の協力・連携を密にし、円滑にプロジェクトを運営していくため、進捗管理等のための会合等を開催すること、2) 外部有識者を招聘した諮問委員会から、本プロジェクトの進捗状況の把握・評価・改善提言・指導等を受けること、を目的としている。また、プロジェクトで得られた成果については、広報や研究活動を通じて積極的に公表するとともに、「富岳」成果創出加速プログラム事業全体のアウトリーチにも協力する。
1)

進捗管理等のための会合に関しては、昨年と同様に効率性を重視した。プロジェクトに参加する
研究者全体の連携を密とすることは必須であるが、運営に過度の労力を割くことは決して得策ではないからである。昨年度と同じく、進捗管理等のための会合を定期的に開催する方針を維持した。これは 2 月に一度の全体会議である（開催日：2021 年 5 月 20 日、7 月 15 日、9 月 17 日、11
月 5 日、2022 年 1 月 31 日、2 月 21 日）。昨年度に引き続き、実施機関が開催し、協力機関を交えて進捗状況や今後の予定等の情報を共有する形で密な連携を図った。

2)

外部有識者を招聘した諮問委員会に関しては、昨年同様、半期に一度開催した。プロジェクト全体の協力・連携を密とし円滑に運営していくことで、より効果的に成果を創出できるよう、諮問委員会は下記の構成としている。諮問委員会からは、本プロジェクトの進捗状況の把握・評価・
改善提言・指導等を受けた。
委員長
平田 直

国立研究開発法人防災科学技術研究所 参与 兼
首都圏レジリエンス研究推進センター長 兼
イノベーション共創本部副本部長

委員
渦岡 良介

京都大学防災研究所 地盤災害研究部門 地盤防災解析研究分野 教授

小野 謙二

九州大学 情報基盤研究開発センター センター長

山中 浩明

東京工業大学 環境・社会理工学院 教授- 110 -

矢来 博司国土地理院 地理地殻活動研究センター 研究管理課長

委員長・委員の専門は、観測地震学・地盤工学・地震工学・数値流体力学・計算科学・物理探査工学・
強震動地震学・固体地球科学である。また、本課題で構築される統合的予測システムを社会実装する際、最大のユーザとなることが想定される国レベルの地震防災・減災の有識者でもある。
第 3 回の諮問委員会は 2021 年 9 月 17 日に Web 形式で開催された。この委員会では、昨年度に指摘のあった本プロジェクトとしての KPI の設定背景や、進捗状況、および今後の進め方について説明した。
KPI に関しては特に、「地形や地下構造、特に沈み込み帯の地震等を扱うには、地下の構造や、地形も入れた解析をすべきという固体地球の分野での世界的な流れがある中で、計算の品質保証が置き去りになっている。この様な現状をふまえ、我々としてはこの様な計算をした結果を、国のハザード対策で活用
いただくことを考えていることから、工学分野でなされているような品質保証の Verification を前提として考え提案したという背景をふくめ、粘弾性応答の計算は品質保証が目的であるために、問題設定として「富岳」全系で計算する必要のないこと、参照解というよりはむしろ収束性等の品質保証など定量的なことが重要である」と考えている点等について説明をおこなった。また、プロジェクトとして、品質保証をしっかり「富岳」を使っておこなったものであることを示していく必要があると考えていることについても言及し、委員から賛同を得られた。成果に係る課題として、「富岳」だからこそできたことをより強調すべきとの指摘を受けた。
第４回の諮問委員会は 2022 年 3 月 3 日に Web 形式で開催された。この委員会では、本年度の研究成果を報告し、評価を受けた。今後に向けて、国際的な位置づけを示すためのエビデンスとして、AGU や EGU
等の国際的な学会での成果発表を通して貢献を示せることが望ましい、企業の人材育成も視野に入れた
具体的なアクション検討を、といった助言を受けた。一方で、企業の若手教育に向けたチュートリアルの取組みや、HPC Asia 2022 における Best Paper 賞と Best Student Poster 賞の受賞については高く評価を受けた。また、最終年度に向けて本プロジェクト終了後の検討を具体化するよう、指摘・助言を受けた。
結論として、本年度の研究の進捗状況・成果に関して、特段の問題がないことが了解された。
プロジェクトで得られた成果については、広報や研究活動を通じて積極的に公表し、今後の展開に資することを進めた。具体的には、国内外の各種学会・研究集会等での発表、および、それら情報を取りまとめた本プロジェクトの和文・英文ホームページ（図 4.1）を介した公表である。ホームページでは見やすさを考慮し、プロジェクトの最新情報をトップページに時系列で掲載するとともに、イベント、研究成果の別で掲載する形式を維持した。情報公開にあたっては、更新が停滞しないよう、可能な限りリアルタイムな情報発信に努め、2022 年 3 月までに昨年度を上回る合計 27 点の記事を掲載した。さらに、
令和３年度成果報告会(2022 年 3 月 3 日 Web 形式)を開催した。
開催形式について、昨年度より引続き Zoom
Webinar を利用すると共に、本年度は YouTube の同時配信を加えた。参加者は全体で 155 名、これは昨年
度の約 1.7 倍増にあたる。成果報告会では本年度の成果に係る報告の後、特別講演として連携機関より代表して 3 名にご登壇いただき、本課題代表者との対談形式も交え HPC への期待や分野を超えた連携等、幅広く闊達な議論が交わされた。

- 111 -

「富岳」成果創出加速プログラム事業全体のアウトリーチへの協力としては、各種学会・研究集会での発表時や、テレビ等のメディア取材時に、「富岳」成果創出加速プログラムおよび本プロジェクトの紹介をおこなった。また、イベントを通した協力としては、2021 年 7 月に開催された「Society5.0 科学博」
による展示、具体的には、本課題の長周期地震動シミュレーションの結果を動画として上映したことや、
2021 年 12 月に開催された「第 14 回 スーパーコンピューティング技術産業応用シンポジウム Society5.0に向けた『富岳』への期待」、2022 年 3 月に開催された「スーパーコンピュータ『富岳』成果創出加速プ
ログラムシンポジウム 富岳百景」での講演等が挙げられる。また、本プロジェクトのホームページ上でもイベント告知等の広報活動をおこなった。
計算機環境については、理化学研究所計算科学研究センターおよび、海洋研究開発機構に整備した計算機環境により、重点課題で開発したアプリケーション群や各関係機関で得られた貴重なデータならびに、本年度の業務で得られた成果等の保管・管理がなされている。今年度は、計算機環境整備の一環として、理化学研究所計算科学研究センターに設置されていた一部機器の海洋研究開発機構

横浜研究所への移設をおこなった。また、研究を推進するための計算機環境の維持・整備として、これら計算機の定常的な監視や定期保守を実施した。

図 4.1.1 和文・英文ホームページイメージ- 112 -

２－３．活動（研究会の活動等）
種別
会議

名称
全体会議（オンライン）

開催日
令和 3 年 5 月 20 日
令和 3 年 7 月 15 日
令和 3 年 9 月 17 日
令和 3 年 11 月 5 日
令和 4 年 1 月 31 日
令和 4 年 2 月 21 日

委員会

諮問委員会（オンライン）

令和 3 年 9 月 17 日
令和 4 年 3 月 3 日

成果報告会

令和 3 年度成果報告会（オンライン）

令和 4 年 3 月 3 日

ワークショップ

チュートリアル令和 3 年 10 月 5 日

２－４．実施体制
業務項目

担当機関

担当責任者

（１）国の被害想定算出のた

神奈川県横浜市金沢区昭和町 3173-25

センター長めの地震・津波災害予測での国立研究開発法人海洋研究開発機構

堀

アプリケーション活用

海域地震火山部門

研究員
縣

高峰
亮一郎

（２）地震に関する災害被害

東京都文京区弥生1-1-1

教授

予測のための大規模アプリケ

東京大学地震研究所

市村 強ーションの開発

計算地球科学研究センター

准教授
藤田 航平

（３）開発した計算コードを

神奈川県横浜市金沢区昭和町3173-25

センター長

実務で活用するための環境整国立研究開発法人海洋研究開発機構

堀備

海域地震火山部門

研究員
縣

高峰
亮一郎

（４）プロジェクトの総合的神奈川県横浜市金沢区昭和町3173-25

センター長

推進

国立研究開発法人海洋研究開発機構

堀

海域地震火山部門- 113 -

高峰

別添１ 学会等発表実績- 114 -