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

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

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

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

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

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

２．令和４年度（報告年度）の実施内容
２－１．当該年度（令和４年度）の事業実施計画（１）国の被害想定のための地震・津波災害予測でのアプリケーション活用
重点課題から開発してきた大規模有限要素コードを用いた長周期地震動の計算結果が、国の被害想定に向けた委員会でこれまで検討対象となってきた。本年度も引き続き、国の被害想定に向けた委員会からの長周期地震動計算結果の公表を目指し、同様な計算での精度検証作業やそのために必要な計算アプリ
ケーションの修正を行う。
（２）地震に関する災害被害予測のための大規模アプリケーションの開発「富岳」の性能を引き出すように、計算科学・計算機科学の最先端技術を駆使することで、地震に関する災害被害予測のためのアプリケーション群の基幹数値解析コンポーネントである、地震・地震動・地殻変動計算等を行う大規模数値解析用コード群に対し、「富岳」実機上でより良好な性能を発揮するように昨年度に引き続き研究開発を行う。本年度は、「富岳」の計算機構の特性を考慮した上で、「富岳」の性能を引き出すように、ソルバーの前処理等の改良を行い、アプリケーション群の性能向上を図る。また、超大規模解析モデルを構築し、これを用いた解析を行うことで、解析手法の有効性を示す。
（３）開発した計算アプリケーション群を実務で活用するための環境整備
国の被害想定算出に使われるものと同等の長周期地震動計算・非線形地盤増幅計算を、建築・土木系の各企業が独自の問題設定で実施するためのサポート体制を、海洋研究開発機構が中心となって整備する。
昨年度、長周期地震動については現実的な地下構造を対象とした実務に近い問題へ適用できるような計算環境の拡充をほぼ完了した。今年度は、プロジェクトに連携機関として入っていただいている企業の方々を全面的に支援して、拡充した計算環境を「富岳」上で用いることによる、実務に有用な問題設定における大規模長周期地震動計算への適用を進める。

-1-

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

うような被害想定においても、差分法による地震動解析を有限要素法によるものに置き換え、幅広い問題設定に対する検討をより精度よく行えるようになることが期待されている。
内閣府政策統括官（防災担当）では、上記のように大規模有限要素法を用いた地震動計算の開発が先行
プロジェクトで進んでいたことに加えて、中央防災会議「東北地方太平洋沖地震を教訓とした地震・津波対策に関する専門調査会報告（平成 23 年 9 月）」において、揺れによる被害を軽減するための対策として、長周期地震動対策等の必要性が指摘されたことや、「南海トラフ沿いの巨大地震による長周期地震動に関する報告（平成 27 年 12 月）」において、震源断層の極近傍に影響地域がある相模トラフは、長周期地震動の高層建築物等の影響評価に至らなかったことなどを政策的な背景として、「相模トラフ沿いの巨大地震に伴う長周期地震動による被害予測の高度化」の課題を進めてきた。その最終段階として、先行プ
ロジェクトの大規模有限要素プログラムをベースとして本プロジェクトで整備した E-wave FEM を用いた
長周期地震動の計算を「富岳」の政策枠で申請し、昨年度より計算を実施してきている。これは、国の被害想定のためのハザード評価の過程に先行プロジェクトからの開発の成果が社会実装されるという点で、本プロジェクトにとって重要な取り組みである。なお、内閣府防災担当は本プロジェクトの連携機関ともなっている。本プロジェクトでは、昨年度に引き続き、内閣府防災担当・業務を担当する防災コンサルタ
ント企業に対し、政策枠を用いた計算を支援した。具体的な支援内容の一部を下記に示す。
● プログラムの改良や使用方法の助言● E-wave FEM の計算精度等を示す計算結果の提供（本年度は②-(c)が該当）● 必要に応じた「富岳」の計算資源やデータ保存領域の拠出● 計算用モデル作成のため、地球シミュレータ（ES4、海洋研究開発機構）の計算資源やデータ保存領
域の拠出、計算実行ファイル等の利用環境整備、及び利用支援
本プロジェクトの支援により実施された長周期地震動計算に基づく被害想定の成果は、別途内閣府防災担当より公表される。
② E-wave FEM の精度検証
本課題においては、①での大規模有限要素法での計算を進めるために必要な精度検証を行い、その結果を①での検討のために提供してきた。なお、本項目に関する数値計算には、「富岳」もしくは海洋研究開
発機構の地球シミュレータを用いた。
E-wave FEM の概要
E-wave FEM で用いられる計算手法について簡潔に述べる。E-wave FEM では弾性波動の支配方程式を有限要素定式化により離散化したシステムを扱う。すなわち
𝐊𝐊𝐊𝐊 + 𝐂𝐂𝐂𝐂 + 𝐌𝐌𝐌𝐌 = 𝐟𝐟である。ここで𝐊𝐊、𝐂𝐂、𝐌𝐌、𝐮𝐮、𝐯𝐯、𝐚𝐚、𝐟𝐟はそれぞれ、剛性マトリクス、ダンピングマトリクス、質量マト
リクス、変位ベクトル、速度ベクトル、加速度ベクトル、外力ベクトルである。ダンピングマトリクスは
レイリー減衰モデル[5]により設定される。レイリー減衰は、指定した周波数の上限値・下限値（FMAX・
FMIN）において、設定したい一定値の Q 値（Quality factor）を近似的に満たすような減衰を導入する。
時間ステップ幅に対する無条件安定性を重視し、陰的な時間積分スキームである Newmark-𝛽𝛽法(𝛽𝛽=1/4,-3-𝛿𝛿=1/2)を適用することで、上記の方程式を𝐀𝐀𝛿𝛿𝐮𝐮𝑛𝑛 = 𝐛𝐛𝑛𝑛 ,

ただし、

𝐀𝐀 =

2
4
𝐌𝐌 + 𝐂𝐂 + 𝐊𝐊,
2
Δ𝑡𝑡
Δ𝑡𝑡

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

𝐪𝐪𝑛𝑛 = 𝐪𝐪𝑛𝑛−1 + 𝐊𝐊𝛿𝛿𝐮𝐮𝑛𝑛 ,𝐮𝐮𝑛𝑛 = 𝐮𝐮𝑛𝑛−1 + 𝛿𝛿𝐮𝐮𝑛𝑛 ,𝐯𝐯𝑛𝑛 = 𝐯𝐯𝑛𝑛−1 +𝐚𝐚𝑛𝑛 = 𝐚𝐚𝑛𝑛−1 −2𝛿𝛿𝐮𝐮 ,
Δ𝑡𝑡 𝑛𝑛4
4𝐯𝐯+
𝛿𝛿𝐮𝐮 .
Δ𝑡𝑡 𝑛𝑛−1 Δ𝑡𝑡 2 𝑛𝑛

毎時間ステップにおいて、巨大な疎行列𝐀𝐀（最も小規模なケースでも 107 オーダーの次元数）を含む線形方程式の求解が必要となる。共役勾配法・マルチグリッド法・精度混合演算を組み合わせた並列 solverを用い、線形方程式の求解を行う[3][4]ことにより、「富岳」や地球シミュレータなどのスーパーコンピ
ュータを使用して多自由度（107-109 のオーダー、あるいはしばしばそれ以上）の問題を高速に解くことが可能となる。
領域端部における反射波は粘性境界条件[6]により処理することを基本としているが、それだけでは処理が十分でない場合がある。そこで、領域端部に吸収帯、あるいはスポンジ領域と呼ばれるバッファ領域
を設け、この領域で弾性波の振幅を人工的に抑制する手法[7]を令和二年度に導入した。この手法では、
(𝑖𝑖)

(𝑖𝑖)

𝐮𝐮𝑛𝑛 = 𝑊𝑊𝐮𝐮0𝑛𝑛

(𝑖𝑖)

(𝑖𝑖)

という形で、弾性波計算の毎ステップにおいて計算量をダンプする。ここで、𝐮𝐮0𝑛𝑛 、𝐮𝐮𝑛𝑛 はそれぞれ、先述した有限要素定式化と Newmark-𝛽𝛽法の組み合わせで定められる更新アルゴリズムによる計算値と、その計算値へのダンピング処理に基づく、第 n ステップでの i 番目の節点での変位である。係数 W は、

-4-

𝑊𝑊 = �exp�−𝛼𝛼(𝐽𝐽0 − 𝑗𝑗)� (𝑗𝑗 ≤ 𝐽𝐽0 )1 (𝑜𝑜𝑜𝑜ℎ𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒)
𝑗𝑗 = 𝑑𝑑 (𝑖𝑖) /𝑑𝑑0

このように定義される。ここで、𝐽𝐽0 , 𝛼𝛼, 𝑑𝑑 (𝑖𝑖) , 𝑑𝑑0 はそれぞれ、吸収帯の厚さに関わるパラメータ、吸収帯内のダンピングの空間変化率を決めるパラメータ、i 番目の節点の最近傍モデル境界からの距離、正規化用距離である。𝑑𝑑0 としては水平方向の最大メッシュサイズを選択する。このようにして、領域端部での反射波が主たる計算領域に影響を与えないようにする。以下の E-wave FEM の計算例では、吸収帯を用いる場合、𝐽𝐽0 ＝120, 𝛼𝛼 = 0.001が使われる。

本課題においては、E-wave FEM について、適切な問題設定下で精度検証を行い、検証の結果に応じてツールの改良も実施してきた。検証の対象としたのは、「波動方程式の求解」、「計算モデルの生成」、「震源断層のモデリング」、「モデル端部における反射波の処理」、「減衰のモデリング」、の 5 項目である。本
年度はこのうち、「波動方程式の求解」について、下記の検討を行った。括弧内は関係する説明のある節を示す。


地形と比較的複雑な地盤構造を伴う構造における長周期地震動計算の差分法との比較（②-a）



地形と現実的な地盤構造を伴う構造における長周期地震動計算の差分法との比較（②-b）



現実的な地盤構造を伴う構造における長周期地震動計算でのメッシュサイズに対する地震動応答変化の確認（②-c）

本課題の実施期間を通じたこれらの精度検証と必要に応じた機能追加により、社会実装に足る地震動計
算を行うための機能と実用性を E-wave FEM が持ち合わせることを確認できたと考えられる。以下に、本年度に実施した検証内容の詳細を述べる。
②-(a) 地形と比較的複雑な地盤構造を伴う構造における長周期地震動計算の差分法との比較
先述のように、国の被害想定のための長周期地震動計算には、差分法による地震動計算手法が用いられてきた。これをE-wave FEMによる計算に置き換えることを検討するにあたり、複雑な地盤構造に基づく問題設定の下で両者を比較することは重要である。地盤構造の複雑さの一要素として、任意形状の地形がある。①で述べた被害想定においては、用いる地盤データの仕様の関係上、地表面を平坦とした計算モデルが扱われているが、今後は現実的な複雑地形を考慮した計算も扱われると見込まれる。有限要素法は、任意形状を非構造要素で表現することができ、さらに自由表面の境界条件を自動的に満たす定式化となっているため、地形を考慮した計算に適している。一方の差分法も、地表面の上に人工的な真空の層を設けることで、地形を近似的に考慮できることが知られている。グリッド幅を細かくすることで、その近似精度をある程度高めることができ、有限要素法とも比較が可能となると考えられる。そこ
で、地形と比較的複雑な地盤構造を伴う人工的な構造において、有限要素法と差分法の数値解を比較し、差異が合理的な範囲にとどまっているかを確認する。ここで人工的な構造を対象とするメリット
は、地盤構造のモデル化以外の両手法の差異が数値解にどう影響するかの切り分けをしやすくなることである。
以下では、有限要素法とその英文名略称であるFEM、差分法とその英文名略称であるFDMを、手法の呼称としてともに用いる場合がある。

-5-

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

 − πf 0 ∆t A(x, y, z ) = exp  Q ( x, y , z ) 
この方法では P 波の減衰 Qp と S 波の減衰 Qs を独立に与えることはできない。また、Q 値の周波数依= Q0 ⋅ f が仮定される。計算領域の側面と下面では、有限の領域を打ち切るために、[10]存性として Qの one-way 波動方程式に基づく無反射境界条件と、[7]の吸収境界条件を併用した境界条件を設けた。

()[7]による吸収境界は、吸収率 W を W = exp − α 2 (J 0 − j ) で表し、境界に近づくに従い指数的に W を増加2させる。α は吸収率の大きさをコントロールする係数、J0 は吸収領域幅を表す。J0 を大きくするほど効
果は大きく、対象とする計算問題に応じて J0 の幅を決める。計算領域の上面では、[11]による
Improved Vacuum Condition を採用し、任意形状の自由地表面を考慮できるようにしている。[12]の方法に基づき、応力の各成分に対して、対応するモーメントテンソルを入力することにより、食い違い震源断層を考慮する。
計算効率の向上のために、領域分割に基づいた並列計算を行った。3 次元の計算領域を南北(NS)、東西(EW)、上下(UD)の 3 方向に分割し、分割した小領域の演算を複数のノードに割り当てて計算負荷を分散させた。計算ステップ毎に MPI を用いた非同期通信を行い、各小領域の袖領域に設けたオーバーラップ
グリッドの速度と応力の値を、隣接するプロセス間で送受信して共有した（例えば、[13][14]）。ノード内の並列計算もプロセス並列とした flat-MPI により行った。
計算に用いた地盤モデルについて説明する。まず、座標系の定義であるが、FDM では+X は北、+Y は東、+Z は鉛直下方をとる。E-wave FEM では+X、+Y、+Z 方向が東、北、上であり、異なることに注意する。地盤モデルの概要を図 1.2.1 に示す（本図は E-wave FEM の座標系に準拠）。このモデルは、[15]に準拠して設定した。出力点として、斜面を跨ぐ Y 方向（図では X 方向）に 1 km 間隔で 14 地点を配置した。また、斜面上および盆地端部の X 方向（図では Y 方向）に、それぞれ 6 地点の出力点を追加で設けた。地盤モデルの物性値を表 1.2.1 に示す。基盤の S 波速度 Vs は 3464 m/s、堆積盆地の Vs は 1000m/s とした。地震波の減衰性を考慮しなかった。
FDM の計算領域は、境界処理に必要な部分を除いて（X, Y, Z）=（17, 17, 10.5）(km) とした。格子間隔は 3 方向とも等しく、昨年度は 20 m、今年度は 10m とした検討を実施した。計算領域端部からの反射波の振幅を小さくするために、端部 250 grid（5 km 幅）で吸収領域を考慮した。時間刻みは0.0005 秒、時間長は 20 秒とした。震源として、（0, 4.5, 4.5）(km)（深さ 3.0km に相当）にダブルカッ-6-プル点震源を与えた。走向、傾斜角、すべり角、地震モーメントはそれぞれ 0 度、90 度、90 度、1018
Nm とした。モーメント・レートの時刻歴関数は、[16]と同様の関数形を用いた。モーメント関数は下記の式で与えられる。
1t−μ𝑀𝑀(𝑡𝑡) = M0 �1 + erf ��� , 𝜎𝜎 = 0.2, 𝜇𝜇 = 4𝜎𝜎2√2𝜎𝜎erf(𝑥𝑥) =2𝑥𝑥� exp(−𝑡𝑡 2 ) 𝑑𝑑𝑑𝑑

√𝜋𝜋 0

M0は地震モーメントである。下の式は厳密には t→∞の極限で 1 となるが、計算の都合上、震源時間関数の計算は t=2s までとした。計算結果に対して、0～2.5Hz までは 1 で、2.5Hz から 3Hz まで線形に変化して 0 となるような周波数特性のローパスフィルターを使い、フィルター処理を施した。

図 1.2.1 に示す 26 地点を観測点とし、計算波形を抽出した。計算には、ES4 を使用し、ベクトルノー
ドで 2048 プロセスの並列計算を行った。計算ステップ数は、40000（＝時間長 20 秒／時間刻み 0.0005秒）であった。グリッド幅 10m の計算には、データの入出力を含めて約 1 時間を要した。
次に、E-wave FEM を用いて同じ比較的複雑な地盤構造における地震動計算を行い、差分法（FDM）を用いた地震動計算結果との比較を行った。先述の通り、E-wave FEM では+X、+Y、+Z 方向が東、北、上である。全体の領域のサイズは FDM の計算とほぼ同様の水平方向約 47 km×47 km、深さ方向 25 km（点震源が Z=21.5 km に）とした。また、領域のサイズは最小メッシュサイズ DS とメッシュ数 NX、NY を指定する形とした。最小メッシュサイズ DS は 50m と 20m の 2 パターンを実施した。Rayleigh 減衰の最大周波数は、ローパスフィルターの形状を考えて FMAX=2.5 Hz、FMIN＝0.25 Hz とした。DS=50m の場
合、総要素数は約 9680 万、総節点数は約 1 億 3200 万となった。吸収帯は X≦12000m, X≧35000m, Y≦
12000m, Y≧35000m, Z≦12000m のいずれかを満たす領域である。DS=20m の場合、総要素数は約 10 億、総節点数は約 13 億 5000 万となった。吸収帯は X≦4800m, X≧42240m, Y≦4800m, Y≧42240m, Z≦4800mのいずれかを満たす領域である。計算時間は 0.01 秒刻みで 15 秒間とした。総ステップ数は 1500 である。所望の震源位置の座標は、このモデル上では（X、Y、Z）＝（28000、23500、21500）(m)となった。
震源パラメータは差分計算と同様である。メッシュ生成および地震動本計算共に ES4 を使用したが、基礎的検討に本課題で供与された Oakforest-PACS の計算資源を用いた。地震動計算については 272MPI プ
ロセスの並列数で実行した。DS=50m では、地震動計算に約 20 分、それ以外の処理に 1 時間あまりを要した。DS=20m では、地震動計算に約 1.7 時間、それ以外の処理に約 2 時間を要した。
FEM（DS=50m）と FDM（グリッド幅 20m）の計算速度波形の比較を図 1.2.2-1.2.4 に、FEM と FDM の波形の差を図 1.2.5-1.2.7 に、Variance Reduction (VR)を図 1.2.8 に示す。すべて FEM の座標系に準拠する。
なお、ここで比較対象とした FDM の結果は減衰無とした場合の結果である。バンドパスフィルターをかける際は、震源時から 15 秒間を切り出し、0.3 秒の cos テーパーを最初と最後に施した。速度フーリエス
ペクトルも同様に 15 秒間の波形を用いた。波形の差、VR およびフーリエスペクトル比については、吸収帯無しの FEM（DS=50m）の結果も示している。また、VR については以下の式により算出した。

-7-

2

∑𝑖𝑖�𝑣𝑣𝐹𝐹𝐹𝐹𝐹𝐹 (𝑡𝑡𝑖𝑖 ) − 𝑣𝑣𝐹𝐹𝐹𝐹𝐹𝐹 (𝑡𝑡𝑖𝑖 )�
� × 100
𝑉𝑉𝑉𝑉 = �1 −∑𝑖𝑖 𝑣𝑣𝐹𝐹𝐹𝐹𝐹𝐹 (𝑡𝑡𝑖𝑖 )2
図 1.2.2-1.2.4 では FEM と FDM の波形はほぼ重なっているように見えるが、差を取りスケールを拡大することによって図 1.2.5-1.2.7 のように差が明瞭になる。FDM と FEM の吸収帯無と吸収帯有との差を比較すると、T=10s 秒後程度以降に目視により差が確認でき、端部からの反射波の影響と考えられる。特に、
図 1.2.1 に示す基盤の斜面上の測線（S11-1 から S11-6 の測線）の UD 成分において反射波の影響による振幅差がそれ以外の部分に比べて相対的に大きい。
VR を見ると、基盤の斜面上の UD 成分において吸収帯を導入することで値がより 100 に近づいている
ことがわかる。斜面上の UD 成分は、元々直達波の振幅が小さく、端部からの反射波の振幅が相対的に大きかったため他の地点・成分に比べ VR が低く出ていたと考えられる。NS 成分については東西測線（S010 から S14-0 の測線）において VR が大きくマイナスとなるため図 1.2.8 には図示していないが、それ以外の観測点の NS 成分の VR は 99 を超えることを確認している。東西測線の NS 成分はラディエーショ
ンパターンからは振幅 0 となる節にあたり、FEM と FDM の震源の入れ方や FDM で食い違い格子を用いていることによる差が強く出ている可能性が考えられる。また、基盤の斜面上の測線（S11-1 から S11-6の測線）の UD 成分の VR が吸収帯を考慮してもなお他の地点の EW・UD 成分の VR より低い理由については、FDM の食い違い格子の影響や斜面の離散化の違いの影響が考えられる。なお、斜面上の測線のNS・UD 成分もラディエーションパターンからは振幅は 0 となる節にあたる。
ここでは DS=50m の FEM とグリッド幅 20m の FDM を比較したが、DS=20m とグリッド幅 10m について比較を行った場合も、質的には同様の議論となる。FEM は非構造要素により地形を精度よく計算に取り入れることが可能である一方、ここで示した人工的な構造においては、FDM においてもグリッド幅を細かく設定することで精度の良い計算が可能となり、両者の結果がよく一致することを確認した。

図 1.2.1 地盤モデルの概要（FEM の座標系に準拠）-8-

表 1.2.1 地盤モデルの物性値
Vp (m/s)Vs (m/s)ρ(kg/m3)

Qp

Qs

堆積盆地

2600

1000

2400

∞

∞

基盤

6000

3464

2700

∞

∞

-9-

盆地内

図 1.2.2

FEM と FDM の計算速度波形の比較（S01-0 から S14-0 の測線）

縦軸は盆地中央（X, Y）＝(23500 m, 23500 m)からの+X 方向の距離（減衰無、FEM 吸収帯有 50m、FDM メッシュ 20m）

図 1.2.3

FEM と FDM の計算速度波形の比較（S04-1 から S04-6 の測線）

縦軸は盆地中央（X, Y）＝(23500 m, 23500 m)からの+Y 方向の距離（減衰無、FEM 吸収帯有 50m、FDM メッシュ 20m）

- 10 -

盆地内

図 1.2.4

FEM と FDM の計算速度波形の比較（S11-1 から S11-6 の測線）

縦軸は盆地中央（X, Y）＝(23500 m, 23500 m)からの+Y 方向の距離（減衰無、FEM 吸収帯有 50m、FDM メッシュ 20m）

盆地内

図 1.2.5

計算速度波形の差（FEM－FDM、S01-0 から S14-0 の測線）

縦軸は盆地中央（X, Y）＝(23500 m, 23500 m)からの+X 方向の距離（黒線：FEM 吸収帯有 50m との差、赤線：FEM 吸収帯無 50m との差）

- 11 -

盆地内

図 1.2.6

計算速度波形の差（FEM－FDM、S04-1 から S04-6 の測線）

縦軸は盆地中央（X, Y）＝(23500 m, 23500 m)からの+Y 方向の距離（黒線：FEM 吸収帯有 50m との差、赤線：FEM 吸収帯無 50m との差）

図 1.2.7

計算速度波形の差（FEM－FDM、S11-1 から S11-6 の測線）

縦軸は盆地中央（X, Y）＝(23500 m, 23500 m)からの+Y 方向の距離（黒線：FEM 吸収帯有 50m との差、赤線：FEM 吸収帯無 50m との差）

- 12 -

図 1.2.8 各測線各成分の VR（黒：FEM 吸収帯有 50m、赤：FEM 吸収帯無 50m）

- 13 -

②-(b) 地形と現実的な地盤構造を伴う構造における長周期地震動計算の差分法との比較
前節で、地形を含む人工的な構造を扱い、有限要素法と差分法の一致度と差の原因について議論した。本節では、地形と実地盤を模した現実的な地盤構造を伴う構造において、有限要素法と差分法の数値解を比較する。昨年度の報告書では、地形なし（平坦）とした現実的な地盤構造における比較を示したが、これに地形を導入したものに相当する。これにより、国の被害想定や実務に今後必要となっていくと考えられるような問題設定下において、両者の差異がどの程度となるかを確認する。
地盤構造データとしては、全国一次地下構造モデル（暫定版）[17][18]の地盤モデルを用いた。海底地形については、日本海洋データセンターJODC による 500 m メッシュ水深データを、陸域地形について
は、国土地理院による 10 m メッシュの DEM データを用いた。標高分布を図 1.2.9 に示す。地盤モデル
は 23 の異なる地震波速度を持つ層からなる。地震基盤（S 波速度 Vs＝3.2 km/s 層）の上面標高分布を
図 1.2.10 に示す。シミュレーションの対象地震は、平成 21 年伊豆半島沖の地震（M5.1）である。図中の●印の位置で、この地震の観測記録が得られている。
水平方向の計算領域は、NS（X）方向 190 km、EW（Y）方向 140 km の範囲（図 1.2.9 および図 1.2.10の図示範囲）とした。深さ（Z）方向では、深さ約 60 km までを計算領域とした。格子間隔は、水平 0.1
km、上下 0.05 km（ただし、深さ 10 km 以深は 0.2 km）とした。計算領域の東西南北の端部および下端に、240 grid（24 km 幅）の吸収領域を設けた。格子点数は（X, Y, Z）=（2380, 1880, 776）grid である（計約 35 億）。計算は、ES4 を用いて、X, Y, Z の 3 方向をそれぞれ 10，8，8 分割した 640 並列により行った。格子間隔が計算結果に与える影響を確認するために、格子間隔を水平 0.05 km、上下 0.02 km（ただし、深さ 10 km 以深は 0.05 km）としたシミュレーションを別途行った。計算領域の東西南北の端部および下端に、480 grid（24 km 幅）幅の吸収領域を設けた。格子点数は（X, Y, Z）=（4760, 3760,。計算は、ES4 を用いて、X, Y, Z の 3 方向をそれぞれ 35，16，8 分割2184）grid である（計約 391 億）した 4480 並列により行った。時間刻みは共通して 0.001 秒、時間長は 200 秒とした。震源として、図1.2.9 の☆印の位置に、図 1.2.11 に示す F-net のメカニズム情報に基づく点震源を与えた。震源時間関数
は、図 1.2.12 に示すライズタイムを 0.9 秒とした二等辺三角形型で与えた。観測記録が得られている地点（図 1.2.9●印）に隣接する格子点で計算結果を出力した。また、海底で計算される地震動特性を調べるために、図 1.2.9▲印の位置に出力点を設けた。計算結果の上下成分は、UP を+、Down を—となるよ
うに符号を変えた。また、計算結果に対して、図 1.2.13 に示す周期特性のバンドパスフィルター処理を行った。
格子間隔を水平 0.1 km、上下 0.05 km とした FDM 計算では、入力データの作成に用いた時間は約 5
時間、プログラム実行時間は約 17 時間、出力データの処理に要した時間は約 0.5 時間である。現状のFDM プログラムでは、将来的に媒質の揺らぎなどを考慮するために、Q 値、格子間隔（NS, EW, UD の3 方向）、媒質番号（1～24）について、全格子点の情報を入力データとして用意し、実行時に、各プロセ
スが割り当てられた計算領域のデータを読み込む形式としている。入力データ量が膨大となるため、入力
データの作成と、プログラム実行時のデータ読込に多くの時間を要する。今後、入力データの削減が必要となる。なお、全国一次地下構造モデル（暫定版）は 23 層から成るが、本検討では地形の考慮のために、空気の層を設けて計 24 層として計算を行っている。
格子間隔を水平 0.05 km、上下 0.02 km とした FDM 計算では、入力データの作成に用いた時間は約- 14 -

60 時間、プログラム実行時間は約 10 時間、出力データの処理に要した時間は約 0.5 時間程度である。この計算では、格子間隔を水平 0.1 km とした FDM プログラムのデータ読込部を修正し、データ読込に要する時間を短縮した。
次に、E-wave FEM を用いて同じ現実的な地盤構造における地震動計算を行い、差分法（FDM）を用いた地震動計算結果との比較を行った。計算領域は、水平方向は図 1.2.9 に示した東西 140 km×南北 190
km の領域に東西南北それぞれに 24 km 幅の吸収帯を加えた東西 188 km×南北 238 km の領域とし、深さ方向は、Z=0 を深さ 108 km とした。深さ方向の吸収帯は Z=0 から Z=72 km までである。観測点数は226 点である。地盤構造の各層境界のメッシュサイズは ds=100m、それ以外の部分は 1 波長当たり 5 要
素以上が生成されるようにメッシュ生成した。FMAX は最小 Vs の 350 m/s および ds から 0.7 Hz とした。全国一次地下構造モデル（暫定版）の物性値のうち Q 値については Qs の値を用いることとした。図
1.2.14 に可視化した第 1 層上端の形状を示す。総要素数は約 28.9 億、総節点数は約 39.0 億となった。
計算はプリプロセスに ES4、地震動本計算およびポストプロセスに「富岳」を用いた。地震動計算は 600
ノード 2,400MPI プロセス数で行った。地震動計算に約 9.5 時間、それ以外の処理に約 18 時間を要した。
E-wave FEM と FDM の比較を行った観測点を図 1.2.15 に示す。FDM の計算は格子間隔の異なる 2 通りの計算を実施しているため、それぞれについて FEM との比較を行った。一般的な格子間隔の場合（水平 100 m 上下 50 m）の FDM（以下、FDM（100m））との、計算速度波形の比較とその差を図 1.2.16、図1.2.17 に、フーリエスペクトルの比較を図 1.2.18 に示す。また、細かい格子間隔の場合（水平 50 m 上
下 20 m）の FDM（以下、FDM（50m））との比較を図 1.2.19～図 1.2.21 に示す。バンドパスフィルターをかける際は、震源時から 200 秒間を切り出し、0.3 秒の cos テーパーを最初と最後に施した。速度フー
リエスペクトルも同様に 200 秒間の波形を用いた。
速度波形（図 1.2.16）を見ると、全観測点を通して、実体波部分と見られる波形の最初の部分は類似しているものの、表面波と見られる後続波の部分において差が比較的大きく、波形の差（図 1.2.17）を取る
ことでそれらが明瞭になる。フーリエスペクトル（図 1.2.18）は、昨年度示した地形なしの際の FEM とFDM（100m）の時に見られたような目立つスペクトルの山谷のずれは見られないが、AAA001 などでは
0.5 Hz に近い部分で山谷が少しずれているのが確認できる。波形の差の傾向や、スペクトル比の傾向は昨
年度の地形なしの場合と変わらないが、フーリエスペクトルは地形なしの場合に比べ同じ FEM と FDM（100m）の差が小さくなっているように見える。地形を考慮に入れたことによりグリッド数が変わり、FEM と FDM のモデルが近くなった可能性が考えられる。FDM（50m）の場合は、差の定性的な傾向は同様であるが、一致度がさらに向上する。グリッド幅を細かくしていくことで、FDM 解が FEM による精度の良い計算結果に近づいていくと考えられる。
これより現実的な構造においても、非構造要素を用いて地形を精度よく取り入れた E-wave FEM による結果とグリッド幅を細かめに設定した FDM の結果がかなりよく一致したといえる。以上の①-(a)(b)のベンチマーキングにより、E-wave FEM と従来法の結果が、一致すべき問題でよく一致し、その差異の程度やそれがどのような部分で出やすいかを調べることができた。

- 15 -

●印：K-NET, KiK-net, MeSO-net 観測点▲印：海底の計算結果確認のために設定した出力点☆印：震央
図 1.2.9 標高分布

- 16 -

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

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

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

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

図 1.2.14 可視化した第 1 層の上面形状（吸収帯は含まない）- 20 -

図 1.2.15 波形ペーストアップに用いた観測点
カラーコンターは、S 波速度 3.2 km/s 層の上面標高を表す。

- 21 -

図 1.2.16 FEM と FDM（100m）の計算速度波形の比較（各波形ペアで振幅を基準化）- 22 -

図 1.2.17 計算速度波形の差（FEM－FDM（100m））（各波形で図 1.2.16 とスケールは共通）- 23 -

図 1.2.18 FEM と FDM（100m）のフーリエスペクトルの比較- 24 -

図 1.2.19 FEM と FDM（50m）の計算速度波形の比較（各波形ペアで振幅を基準化）- 25 -

図 1.2.20 計算速度波形の差（FEM－FDM（50m））（各波形で図 1.2.19 とスケールは共通）- 26 -

図 1.2.21 FEM と FDM（50m）のフーリエスペクトルの比較- 27 -

②-(c) 地形を含む現実的な地盤構造における長周期地震動計算でのメッシュサイズに対する地震動応
答変化の確認有限要素法や差分法などの数値計算においては、計算用のメッシュサイズやグリッド間隔等の離散化
幅を小さくするほど計算精度が向上するが、それに従って必要な計算コストも大きくなる。そのため、適切な離散化幅を選択することが重要となる。長周期地震動の差分法計算においては、対象とする周期の下限における 1 波長を 6～10 グリッド（E-wave FEM で用いられる四面体二次要素では 3～5 要素に相当）以上で離散化できるように設定されることが多い。しかし実際には、ある精度を担保するために必要な離散化幅は問題設定依存であることが指摘されており、本課題における昨年度までの検討でも同様のことを確認している。また、数値計算に由来する誤差が十分小さいかどうかの絶対的な判断は不可能であり、計算の用途に照らし合わせることで初めて意味のある議論ができるようになる。他方、国の被害想定においては、E-wave FEM を用いた長周期地震動計算の際、「1 波長当たり 5 要素確保」という条件を基準として計算モデルの作成がなされてきた。先述のような理由から、本設定が用途に対して十分な計算精度を確保するものとなっているかを確認することが望ましい。そこで、国の被害想定で用いられるものと同等の長周期地震動計算設定における E-wave FEM の計算結果に対し、非構造格子のメッシュサイズに対する数値
解の変化を確認する。着目するのは、結果公表の対象となる擬似速度応答スペクトルである。その変化の程度を、国の被害想定で考慮される震源パラメータの多様性に基づく多数ケース計算における擬似速度
応答スペクトルのばらつき度合いと比較し、議論する。なお、ここで議論する計算結果は、想定される相模トラフ巨大地震に対する関東地方の地震動強さという重要な情報を含むため、応答の絶対値を示すことは避け、異なる設定に対する結果どうしの相対的な差についてのみ検討する。
精度確認の対象とするE-wave FEMの計算設定においては、関東平野を含む東西方向182km、南北方向
168km、鉛直方向140kmの計算領域に、4つのアスペリティからなる震源断層モデルが配置されている（図1.2.22）。この震源モデルは相模トラフのプレート境界面に位置しており、いわゆる「相模トラフ
巨大地震」を模したものである。最西部に位置するアスペリティから破壊が始まり、東方にあるアスペ
リティに徐々に伝播していく。この震源断層モデルは、多数ケースの震源パラメータを考慮したうちの1ケースに相当する。地震波速度構造は、防災科学技術研究所のJ-SHISデータベース[19]をもとに、地震波速度が一様な41の層からなる三次元不均質速度構造として設定されている。J-SHISの地盤データに対してこれまで行われてきたValidationの内容を踏まえて、地形は平坦としている。今回検討対象とす
る国の被害想定と同等の基準設定ケース(以下BASEと呼ぶ)においては、非構造要素メッシュを1波長当たり最低5要素確保するように生成する。この計算用モデルを用い、周期2秒以上の帯域に対して、300
秒分（時間刻み幅0.02秒×15,000ステップ）の地震動が計算されている。事前の解析により、この時間
刻み幅が十分小さいことはすでに確認済みである。比較対象として、1波長当たり最低10要素確保する
メッシュを準備した。このケースをNOW10（The number of wave is 10）と呼ぶ。NOW10のモデル規模は比較的大きくなり、約156億自由度、約39億要素となった。その他の設定はBASEと同じものとし、地震動を計算した。NOW10を用いた計算には、「富岳」の1,024計算ノードを用いて約7時間を要した。図
1.2.23に、計算用モデルの概観と計算された地震動のイメージ図を示す。
BASE と NOW10 の地震動計算結果を、擬似速度応答スペクトルについて比較する。応答スペクトルは、ある固有周期 T をもつ一質点系モデルが地震動を受けた際に生じる応答（加速度、速度、変位）の最大値- 28 -

である。擬似速度応答スペクトルは加速度応答スペクトルの最大値に T/2πをかけることで求められる。
国の被害想定においては、長周期地震動の評価結果として、水平 2 成分の時刻歴応答合成をもとにした
擬似速度応答スペクトルが用いられるため、ここでもそれに倣うこととした。下記のような手順により、両ケースにおける擬似速度応答スペクトルを網羅的に比較した。なお、応答スペクトルの計算には、[20]のプログラムを用いた。
1) 計算対象となる地表面全体を覆うように、990 点の擬似速度応答スペクトル評価点（図 1.2.22 のマ
ゼンタ点）を格子状に配置する。
2) BASE と NOW10 双方の地震動計算結果から、990 の評価点における擬似速度応答スペクトルを、周期
2、4、6、8、10 秒について計算する。
3) 固有周期 5 ケースごとに、990 点のすべてについて、擬似速度応答スペクトル(pSv)の相対差を計算する。相対差は、(pSvBASE-pSvNOW10)/pSvNOW10 と定義する。
4) 固有周期 5 ケースごとに、擬似速度応答スペクトル相対差に対するデータ数 990 のヒストグラムを作成する。
作成されたヒストグラムを図 1.2.24 に示す。相対差が 0％をピークとして対称的に分布しており、ほとんど点で±5％より小さいレベルに収まっていることがわかる。周期 2、4 秒の場合には、±5％を超える点が数十点見られるが、±10％を超えるのは数個のみである。国の被害想定では、震源パラメータの多様性を考慮し、破壊開始点や破壊伝播速度に揺らぎを持たせた多数ケース計算の統計量を検討に用いる。
この多数ケース計算における擬似速度応答スペクトルの相対的なばらつきの範囲は 20～30％程度となっている。BASE と NOW10 での擬似速度応答スペクトルの相対差はほとんどの評価点でこれより十分小さい
レベルとなっていることが確認された。NOW10 とメッシュをさらに細かくしたモデルでの計算結果の差はこれより大幅に小さくなると考えられるため、「BASE での 1 波長当たり 5 要素確保」という設定が、国の被害想定の用途に対しては十分であることを確認した。
ここで、どのような条件で差が大きくなっていたかを考察する。相対差が 5％以上と比較的大きくなる
評価点は、周期 2、4 秒の短周期の場合に多く含まれている。このことは BASE と NOW10 の短周期の帯域において計算時刻歴波形の差が比較的大きくなっていることを示唆する。短周期成分の計算により細かいメッシュが必要となることを考えると、この結果は妥当と考えられる。次に、周期 4 秒を例に、評価点の場所による相対差の違いを見てみると、差の大きい地点は千葉県北部から茨城県南部付近に比較的集まっているように見える（図 1.2.25）が、計算モデルの空間変化やその BASE と NOW10 での違いとの対応
関係は、このままでは明白ではない。そこで、相対差が大きかった観測点の一つと破壊開始点を結ぶ測線に注目する。測線の位置を図 1.2.25 に示す。この線は、周期 4 秒における擬似速度応答スペクトルの相対差が 9.5％程度と大きい C という点を通っている。比較のため、同じ線の付近にあり比較的相対差の小さい点 A、相対差が点 A と点 C の中間程度である点 B と合わせ、各地点での計算加速度波形を図 1.2.26に示す。どの点の波形も全体としては BASE と NOW10 はよく一致しているように見える。ただ、点 C では、
170 秒くらいから、後続波での振幅・位相双方のずれが目立ってくる。擬似速度応答スペクトルの比較図を図 1.2.27 に示す。点 C の周期 4 秒付近における BASE と NOW10 の差が、他の点・帯域と比べるとかな
り大きくなっていることがわかる。これらの差と計算モデルの関連を調べるため、モデル化された P 波
速度構造の測線下での断面を BASE と NOW10 で比較してみる。全体図を図 1.2.28 に、点 A、B、C 付近の拡大図を図 1.2.29 に示す。全体図では両者の差は視認できないが、拡大図においては、両計算モデルに- 29 -

おける浅部の地盤構造のモデルの差が明らかである（縦方向を 50 倍程度伸長しており、誇張された差が見えていることに注意）。差の小さい点 A より震源側（図左方向）の地下においては、地盤構造が水平成層に近い単純な構造となっており、メッシュサイズを変更したことが及ぼす計算モデルへの影響が小さ
いことがわかる。一方、差がやや大きい点 C と中間的な点 B の地下においては、深さ 200m 程度までの地盤構造が複雑となっており、メッシュサイズの計算モデルへの影響が大きくなっている。これらは周期 2、4 秒といった短周期側の地震波に影響を与えるスケールの構造なのだと考えられる。このように、局所的に BASE と NOW10 の擬似速度応答スペクトルの差が大きくなる場合があることの主要因の一つが、浅部の地盤構造が複雑な地域においてメッシュサイズを変更したことが及ぼす計算モデルの違いである可能性が高いことが、計算モデル断面図の可視化・比較から示された。

図 1.2.22 関東地方計算領域図
白線、黒点線、マゼンタ点、色コンターはそれぞれ、海岸線、震源モデル位置、計算
値評価点、地震基盤面深さを示す。

- 30 -

図 1.2.23 NOW10 での地震動計算
上：対象とする関東地方に対して生成した有限要素モデルの概観。色コンターは地震基盤の深度に対応しており地震基盤以浅の層は透過処理してある。白線は海岸線の位置を示す。
下：計算された波動場の可視化結果のスナップショット。深度に対するグレースケールコンター付き地
震基盤形状の上に、地震基盤以浅の層が半透明で示される。色コンターは z 方向変位で、赤が正の値、青が負の値を示す。白線は海岸線の位置を示す。右上の数字は破壊開始からの経過時間(分：秒)。

- 31 -

図 1.2.24 各固有周期での水平合成擬似速度応答スペクトル(pSv)の、BASE と NOW10 の相対差(%)のヒストグラム(サンプル数 990、階級数 48)
上行左から、周期 2，4，6 秒、下行左から周期 8, 10 秒のもの。

- 32 -

図 1.2.25 周期 4 秒における擬似速度応答スペクトルの相対差（絶対値）の空間分布
破壊開始点付近から相対差の大きい点 C を結ぶ測線（黒鎖線）と、その付近にある点 A、B の場所（アルファベットはその左下の点の名称）も示す。灰鎖線は震源モデル位置。

- 33 -

図 1.2.26 測線上の評価点における計算加速度波形（遮断周波数 0.5Hz のローパスフィルターで処理済。点ごとに応答の最大値で正規化済）
上から、点 A、B、C での x、y、z 方向の波形。黒線と赤鎖線はそれぞれ BASE、NOW10 を示す。

- 34 -

図 1.2.27 測線上の評価点における水平合成擬似速度応答スペクトル（点ごとに応答の最大値で正規化済）注目する 5 つの周期だけでなく、周期 2～20 秒までの帯域で、周期 0.2～0.5 秒ごと（刻み幅は周期帯によって変化）のスペクトル値を示している。上から、点 A、B、C での値。黒線と赤鎖線はそれぞれ
BASE、NOW10 を示す。
- 35 -

図 1.2.28 図 1.2.25 で示す測線下の断面における計算モデルでの P 波速度構造
上が BASE、下が NOW10 の計算モデルを示す。赤枠は図 1.2.29 における拡大部分。
縦軸は計算モデル上の z 座標であり、140,000m が地表面に相当する。

- 36 -

図 1.2.29 図 1.2.28 で示される P 波速度構造の赤枠部分の拡大図と点 A、B、C の位置
上が BASE、下が NOW10 の計算モデルを示す。縦軸は計算モデル上の z 座標であり、140,000m が地表面に相当する。縦方向を 50 倍程度伸長しており、誇張された差が見えていることに注意。

- 37 -

参考文献[1] Hisada, Y., and S. Yamamoto. "One-, two-, and three-dimensional site effects in sedimentfilled 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

groundmotionmodellinginbasinsincluding

topography 、Geophysical Journal International 、 Volume 177 、 Issue 3 、 Pages 1221 1232 、June 2009[3] Ichimura, Tsuyoshi, et al. "Physics-based urban earthquake simulation enhanced by 10.7
BlnDOF× 30 K time-step unstructured FE non-linear seismic wave simulation." SC'14. IEEE,2014.
[4] 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 theInternational Conference for High Performance Computing, Networking, Storage and Analysis(pp. 15-26). 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] Lysmer, J., & Kuhlemeyer, R. L. (1969). Finite dynamic model for infinite media. Journalof the Engineering Mechanics Division, 95(4), 859-878.
[7] Cerjan, C., Kosloff, D., Kosloff, R., & Reshef, M. (1985). Discrete acoustic and elastic
wave equations. Geophysics, 50(4), 705–708.
https://doi.org/http://dx.doi.org/10.1190/1.1441945[8] 林田智宏 , 岡元太郎 , 竹中博士 : 速度 応力型スタ ガード格子差分法を用いた 2 次元及び 3 次
元地震波動計算コードの作成 , 九州大学理学部研究報告 地球惑星科学 , 20 巻 , 3 号 , pp. 99 110,
1999.
[9] Graves, R.: Simulating seismic wave propagation in 3D elastic media using staggered-gridfinite differences, Bull. Seism. Soc. Am., Vol. 86, No. 4, pp. 1091-1106, 1996.
[10] 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.
[11] Zeng, C., Xia, J., Miller, R., Tsoflias, G.: An improved vacuum formulation for 2D finitedifference

modelingofRayleighwaves

includingsurface

topography

andinternal

discontinuities, Geophysics, Vol. 77, No. 1, pp. T1-T9, 2012.
[12] Pitarka A.: 3D Elastic finite-difference modeling of seismic motion using staggered gridswith nonuniform spacing, Bull. Seismol. Soc. Am., Vol. 89, No. 1, pp. 54-68, 1999.
[13] 笠松健太郎, 山中浩明: 中小地震の地震動シミュレーションによる 2004 年新潟県中越地震震源域の深部地盤構造モデルの構築, 物理探査, 59 巻, 5 号, pp. 475-484, 2006.

- 38 -

[14] 笠松健太郎, 加藤研一: 大規模並列計算による 3.11 本震の長周期地震動シミュレーション, 地球惑星科学連合大会, 2016.
[15] 小林広明，渡辺哲史，笠松健太郎，加藤研一，縣亮一郎，堀高峰：富岳を用いた大規模地震動計算に向けた Verification－地形と盆地を伴う構造でのベンチマークテスト―，日本地球惑星科学連合 2022
年大会，STT41-01，2022.
[16] 吉村智昭, 永野正行, 久田嘉章, 青井真, 早川崇, Citak, S., 松島信一, 川辺秀憲, 上林宏敏: 強震動
予測手法に関するベンチマークテスト：数値解析手法の場合（その 2）, 18 巻, 38 号, pp. 95-100,
2012.
[17] Kazuki Koketsu, Hiroe Miyake, Yasuhisa Tanaka, et al. A proposal for a standard procedureof modeling 3-D velocity structures and its application to the Tokyo metropolitan area,
Japan. Tectonophysics, 472(1- 4):290–300, 2009.
[18] 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.
[19] Japan Seismic Hazard Information Station (J-SHIS) [Online]. http://www.j-shis.bosai.go.jp/.
[20] 大崎順彦：新・地震動のスペクトル解析入門，鹿島出版会，1994（２）地震に関する災害被害予測のための大規模アプリケーションの開発
本項では、計算科学・計算機科学の最先端技術を創出しこれを駆使することで、地震に関する災害被害予測のためのアプリケーション群の基幹数値解析コンポーネントである、地震・地震動・地殻変動計算等を行う大規模数値解析用コード群に対して、スーパーコンピュータ「富岳」実機上でより良好な性能を発揮するよう開発した研究について述べる。
本年度は、「富岳」の計算機構の特性を考慮した上で、その性能を引き出すように、ソルバーの前処理等の改良を行い、アプリケーション群の性能向上を図った。また、超大規模解析モデルを構築し、これを用いた解析を行うことで、解析手法の有効性を示した。具体的には、「①地殻モデルの不確実性を考慮した超大規模地殻変動解析」、「②シミュレーション内で生成されるデータの学習による高詳細粘弾性地殻変動解析の高速化」、「③「富岳」で開発した高効率行列ベクトル積手法の汎用化のための開発」、「④地震観測データをシミュレーションにて有効活用するための高速データ解析手法の開発」を行った。なお、①は SC22 にて Gordon Bell Prize Finalist に選ばれるなど高く評価されている。また、②、③は SC22 のワークショップに採択され、④は International Conference on Computational
Science において extended paper（発表された査読付き論文のうち優れた論文が選定される）に選ばれるなど、それぞれ高い評価を得ている。
以下に、上記の①～④の具体的な成果をまとめる。
① 地殻モデルの不確実性を考慮した超大規模地殻変動解析①-(a). はじめに
自然災害のひとつである地震は、世界中の様々なところで生じる[1]。大地震は、人命・財産・社会に対- 39 -

し時に破滅を与える、人類にとって大きな脅威である。実際、例えば、平成 23 年の東日本大震災では、
2 万人以上の死亡・行方不明者が発生し、平成 27 年時点でも 20 万人以上が避難生活を余儀なくされた。
そしてそれ以上の被害を及ぼす懸念があるものとして、例えば、日本の南海地方で発生が懸念される南海トラフ地震があり、また、同様の形態の地震としてアメリカ西海岸の Cascadia Subduction Earthquake
やニュージーランド北島東海岸の Hikurangi Subduction Earthquake などがある。
大地震のプロセスは、
i)プレートの動きによりプレートにひずみがたまり限界が来ると断層を生じ、ii)断層により生じた地震動
が地殻内を伝わり、iii)さらに表層近傍の軟弱層で強められた地震動が人命・財産・社会に被害を及ぼす(図 2.1.参照)。また、i)は大規模な津波を引き起こす場合もある。地震被害軽減は人類の宿願であり、このような社会的要請が強い地震問題へのスーパーコンピューティング分野からの貢献が強く期待されてきており(e.g. [2])、様々な成果が創出されている。例えば、ii)のための手法として[3], [4], [5], [6]が、iii)のための手法として[7], [8], [9]などが Supercomputing Conference (SC)で発表されており、地震問題に対する解析能力の向上に大きな貢献を果たしており、国や地域の防災能力向上への貢献が期待されている。

図 2.1. 地震プロセスは、i)プレートの運動・地震発生、ii)地震動伝播、iii)地盤増幅・都市及び社会の応答までをも含む広大な領域を対象とする。なお、i)プレートの運動・地震発生は、場合によっては大規模
な津波も引きおこす。これらのシミュレーションの領域は広くかつ分解能は非常に細かいため、解析コストは膨大となる。その中でも、本研究で対象とする物理ベースシミュレーションによる i)プレート運動・
地震発生の解析コストは膨大となる。

一方、i)に関するスーパーコンピューティングにおける検討はまだ少ない。i)は地震発生の時期・場所・規模を予測する地震予測や、どのような地震・津波が起きうるのかという震源シナリオを構築する上で重要な役割を持ち、ii)及び iii)の入力となる。そのため、地震被害軽減を目指す ii)iii)の貢献の大きいことは論をまたないが、i)の信頼性向上はそれを凌駕する可能性がある。例えば、近年の海陸での観測網の発達により、大地震前のゆっくりした断層すべりなどが捉えられたという報告がなされている（例えば[10]など）。ゆっくりすべりが発生したからといって必ずしも大地震になるとは限らないが[11]、その可能性が普段よりも高いとは言える。ゆっくりすべりによる近傍の大地震発生確率の変化を物理モデルに基づいて推定する試みもニュージーランドで行われており[12]、統計モデルによる確率評価と並んで政策判断にも活用された[13]。こうした物理モデルにもとづいた解析を行う際、従来はシンプルな幾何形状のモデル- 40 -

が用いられていたが、最近は海底地形や沈み込むプレートの形状などの３次元的な地下構造をモデル化した解析に移行しつつある（例えば[14]など）。ここで問題になるのが、複雑な地下構造情報の曖昧さが解析結果に与える影響である。その影響を定量的に評価するためには、複雑で三次元的な幾何形状を持つ
地殻構造の情報の曖昧さを考慮しつつ、十分広い領域・十分高い分解能での地殻変動解析により i)を解くことが必要となる。しかし、そのためには超大自由度の方程式を解く必要があり、SC のゴードンベルセッションにおいてもまだ試みられていない。本研究では、i)のベースとなる地殻構造の情報の曖昧さを考慮した地殻変動解析を実行可能な、32 兆自由度までの超大規模確率有限要素解析を「富岳」全系スケー
ル で 実 現 す る ソ ル バ ー ACROSS を 開 発 し た 。 従 来 法 で は 難 し い と さ れ て き た Uncertainty
Quantiﬁcation を、高詳細地殻変動解析で考慮することが出来るようになるため、ACROSS は地震被害軽減へ向けた大きな成果となる可能性が期待される。また、昨今注目されている Slow Earthquake [15]等地震の物理を理解する上でも Uncertainty Quantiﬁcation は必須であることから、サイエンスへの貢献も期待される。なお、従来の物理空間に加え、確率空間をも同時に離散化の対象とする本研究のようなアプローチは超大規模自由度となりすぎるため、他のどの分野をみても例は少ない。一方で、本研究と同様に Uncertainty Quantiﬁcation は多くの分野で必要とされていることから、このような方向性の研究は今後 Capability Computing の一分野として、大規模シミュレーションの貢献が望まれる新しいフロン
ティアになる可能性がある。
①-(b). 関連する研究のなかでの本研究の位置づけ
i) の地殻変動解析は、線形/非線形・静的/粘弾性でモデル化される。本研究では、線形/非線形・静的/
粘弾性解析のベースとなる静的線形弾性を対象として開発を進めた（本手法を拡張することで、線形/非線形・静的/粘弾性を考慮した解析は実現可能である）。本問題は、地表での応力フリーの境界条件を解析的に満足し、地殻構造の三次元的で複雑な形状と物性の高い heterogeneity を適切に取り扱う必要がある
ことから、105−6 × 105−6 × 105−6 m オーダーの領域𝑉𝑉で、
�𝑐𝑐𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 𝑢𝑢𝑘𝑘,𝑙𝑙 (𝒙𝒙)� = 𝑓𝑓𝑗𝑗 (𝒙𝒙)
,𝑖𝑖

in

𝑉𝑉

(1)

を変位関数𝑢𝑢𝑖𝑖 について解く低次有限要素による有限要素解析で解くことが望まれる。なお、側面と底面へ境界条件として半無限境界を設定することから、通常の固定境界条件等と比較して収束性は悪化する。ここで、𝑐𝑐𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 , 𝒙𝒙, (,𝑖𝑖 ), 𝑓𝑓𝑖𝑖 は、弾性テンソル、空間座標、𝑖𝑖方向への偏微分、外力である。解の収束性を担保するために必要な離散化分解能からすると、この方程式の未知数の自由度は 1012−15 程度となるため、結果的に、implicit iterative solver による大自由度倍精度非構造四面体二次要素有限要素解析を開発することとなる。一方で、昨今の観測データの充実を踏まえても地殻構造の情報には曖昧さが残っており、シミュ
レーションの品質保証のためには、これを考慮した Uncertainty Quantiﬁcation が重要となるが、これだけの規模の解析では難しいとされてきた。
本研究では、以降の節で詳細を説明する超大規模確率有限要素解析を実現することにより、地殻変動解
析の Uncertainty Quantiﬁcation を可能とした。本研究により実現された SoA（state-of-the-art）として、従来は不可能とされていた、前述の南海トラフ地震での検討例を一つ示す。巨大地震が発生し社会が混乱に陥っている時、地下ではそれに引き続いて、同等かそれ以上の地震に向けた準備が進んでいることがある。大地震が発生するとその周囲にゆっくりした断層すべりが広がり（afterslip）、それが広がった- 41 -

先で応力集中を起こし、壊れ残っていた範囲での地震（同程度、あるいはより大きな規模）を誘発するのである[16]。事実、南海トラフ地震がその典型例であり、マグニチュード 8 クラスの地震の１日半後ある
いは２年後に同程度以上の地震が隣接領域で発生している。このように、巨大地震が発生した後、連鎖的に同程度あるいはより大きな地震が発生する危険性を把握するために、巨大地震発生とその後の推移を
曖昧さ込みでモニタリングし、引き続く巨大地震の発生可能性を評価することの重要性は、社会の安全を守る観点から非常に高い。本研究における適用例では、実際に大地震が過去に誘発されたことのある南海
トラフ地震の例を模擬する。ここでは図 2.2.a）に示す領域について、大地震の続発を含めた過去の地震
発生の繰り返しパターンを再現したシミュレーションを行った（地震時すべりおよび afterslip [17]に対する南海トラフ領域の Uncertainty Quantiﬁcation を含んだ地殻変動の推移予測を行った）。既往研究（[18], [19], [20]）での各層が取りうる値の上限・下限により、本解析においては各層の体積弾性率 K とせん断剛性率 G について曖昧さを導入している(図 2.3.)。本解析では、全域において、ポイントワイズで地殻変動の確率応答分布が求められている。このすべての点で得られた確率応答を処理することで、地震直後及び 7 日後の地殻変動結果の 3σ、平均、標準偏差を求めた（図 2.2.c））。物性の曖昧さにより応答の振幅が大きく異なっていることがわかり、Uncertainty Quantiﬁcation を行うことの重要性がわかる。地震時に大きな変位が出ているところでは平均と標準偏差の比率はそこまで大きくないが、地震時に大きな変位がみられるエリアの周辺部、また、afterslip の場合には、相対的に平均と標準偏差の比率が大きくなっている。隣接する領域での大地震の誘発の評価にとって鍵となるのは、すでに起きた地震時のすべりの周辺部や afterslip のすべり量やその広がりによって生じる応力変化である。すでに、これらの領域の真上である水深 1900-4300 m の海底に先端的な地殻変動観測システム[21]（観測点分布を図 2.2.a）に示す）が展開されており、本適用例結果とあわせることで、大地震が発生した場合にこれらの領域のすべ
りや応力変化を、より適切に推定することが期待されている。
Uncertainty Quantiﬁcation が可能な方法の典型的な手法は確率有限要素法であるが、現時点では、提案手法である ACROSS 以外上記の適用例規模を求解可能なものはなく、比較可能な SoA はない。代替的に、上記の適用例を実現可能な方法として、deterministic な非構造四面体二次要素有限要素解析を用いたモンテカルロシミュレーションが考えられるので、これとの比較を行う。適用例を deterministic に計算した場合の有限要素モデルの自由度は 111,196,206,711 (水平要素サイズ 375 m)となる。この規模の地殻変動解析を実行可能な SoA は[22]を拡張し Summit 用に開発した[23,24]と考えられるが、直方領域の等方弾性のみ、かつ、一部の GPU でしか有効でない方法なので、SoA として適当でない。そこで、[23,24]で複雑形状・等方弾性体に限らない一般の弾性体を対象とした汎用アルゴリズムについて
Summit 用にチューニングした手法を用いてこの有限要素モデルを用いたモンテカルロシミュレーショ
ンを行った場合の解析コストが SoA となる。適用例を模した次節の検証問題程度の試行回数（104）が少なくとも必要だとして、性能比較情報がある全系性能計測問題同士で比較すると、Summit 全系 4544
node を活用して 1300 億自由度の問題を 104 回解くには 52451 s かかることとなる。一方で、ACROSS
では、「富岳」全系 147456 ノード利用時で同等の問題が 234 s で求解できている。つまり、ACROSS は、
Summit 全系で SoA を実行した際と比べて 224 倍の高速化を達成していることが分かる。なお、この結
果は massive capability computing が capacity computing よりも桁違いに高速になりうるという一例であり、超大規模計算機環境の整備の有効性が示される計算例となっている。

- 42 -

図 2.2. 「富岳」全系による、領域サイズ 2496×2496 km、最小要素サイズ ds = 375 m、
3𝑛𝑛 × 𝑛𝑛𝑃𝑃𝑃𝑃 = 32.0兆自由度の Uncertainty Quantification を含んだ地殻変動解析例

図 2.3. 南海トラフ地震の解析に使用した有限要素モデルと物性- 43 -

図 2.4. モンテカルロシミュレーションと提案手法により得られる確率応答分布の比較
a）有限要素モデル((16.5, 20.5, 17.5) km 地点で点震源を入力)、b) 材料特性、c) (18500, 11000, 24956.2) m での変位量の比較。

①-(c). 開発手法
本研究では、地殻構造を成層構造として𝐾𝐾と𝐺𝐺に各々曖昧さがある地殻変動解析を対象としているが、定式化が煩雑であるため、簡単のため、ヤング率に曖昧さがある場合の確率有限要素法（e.g., [25]）の定式化を示す。なお、これは𝐾𝐾と𝐺𝐺に各々曖昧さがある定式化と本質的に同等である。地殻をヤング率𝐸𝐸が確率変数である線形連続弾性体𝑉𝑉とする。𝑉𝑉の側面と底面には無限境界が適用され、外力𝑓𝑓𝑖𝑖 (𝒙𝒙)が与えられた
場合、変位𝑢𝑢𝑖𝑖 は次の確率境界値問題の解となる。

�𝐸𝐸(𝒙𝒙, 𝜔𝜔)𝑐𝑐̅𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 𝑢𝑢𝑘𝑘,𝑙𝑙 (𝒙𝒙, 𝜔𝜔)� = 𝑓𝑓𝑗𝑗 (𝒙𝒙) in 𝑉𝑉 × Ω
,𝑖𝑖

(2)

ここで𝜔𝜔は標本空間Ωの標本点を表し、𝒙𝒙は連続体𝑉𝑉内における点の位置座標を表している。𝑐𝑐̅𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 は、ポア𝜈𝜈ソン比𝜈𝜈によって𝑐𝑐̅𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 = (1+𝜈𝜈)(1−2𝜈𝜈) 𝛿𝛿𝑖𝑖𝑖𝑖 𝛿𝛿𝑘𝑘𝑘𝑘 +12(1+𝜈𝜈)(𝛿𝛿𝑖𝑖𝑖𝑖 𝛿𝛿𝑗𝑗𝑗𝑗 + 𝛿𝛿𝑖𝑖𝑖𝑖 𝛿𝛿𝑗𝑗𝑗𝑗 )として与えられる。

物理空間と確率空間で𝐸𝐸と𝑢𝑢𝑖𝑖 を離散化する確率有限要素法により Eq. (2)を解く。𝑉𝑉の一部である𝑉𝑉𝑘𝑘 内で���𝑘𝑘 かつ標準偏差𝜎𝜎𝑘𝑘 の曖昧さを持つ場合を考え、互いに独立な標準正規分布（期待値 0、分散 1
𝐸𝐸が均一値𝐸𝐸の正規分布）に従う確率変数である標準正規確率変数𝜉𝜉𝑘𝑘 を用いて、
𝑛𝑛

𝐾𝐾𝐾𝐾 (𝑘𝑘)𝐸𝐸 (𝒙𝒙)𝜉𝜉𝑘𝑘 ,𝐸𝐸 = 𝐸𝐸 (0) (𝒙𝒙) + ∑𝑖𝑖=1とする。𝑛𝑛𝐾𝐾𝐾𝐾 は𝑉𝑉𝑘𝑘 の数、𝐸𝐸 (0) は点𝒙𝒙における平均であり、- 44 -(3)

(4)𝐸𝐸 (𝑘𝑘) (𝒙𝒙) = 𝜎𝜎𝑘𝑘 𝐻𝐻𝑘𝑘 (𝒙𝒙) (𝑘𝑘 = 1,2, … , 𝑛𝑛𝐾𝐾𝐾𝐾 )である。𝐻𝐻𝑘𝑘 (𝒙𝒙)は𝑉𝑉𝑘𝑘 内で値 1、その他の領域で値 0 をとる Heaviside 関数である。次に、 𝑢𝑢𝑖𝑖 は物理・確率空間で連続的に分布すると期待されることから、Polynomial Chaos (PC)展開により以下のように展開す
る。
𝑛𝑛

(𝑝𝑝)

(5)

𝑃𝑃𝑃𝑃
𝑢𝑢𝑖𝑖 Ψ𝑝𝑝
𝑢𝑢𝑖𝑖 = ∑𝑝𝑝=0(𝑝𝑝)

ここで、Ψ𝑝𝑝 は PC 展開の基底、𝑢𝑢𝑖𝑖 はその係数である。

最終的に、確率境界値問題 Eq. (2)から導かれる確率変分問題である汎関数、
1∫𝑉𝑉×Ω 2 𝐸𝐸(𝑥𝑥, 𝜔𝜔)𝑐𝑐𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 𝑢𝑢𝑖𝑖,𝑗𝑗 (𝑥𝑥, 𝜔𝜔)𝑢𝑢𝑘𝑘,𝑙𝑙 (𝒙𝒙, 𝜔𝜔) − 𝑢𝑢𝑖𝑖 𝑓𝑓𝑖𝑖 𝑑𝑑𝒙𝒙𝑃𝑃(𝑑𝑑𝑑𝑑),(6)に対して、Eq. (3)及び Eq. (5)を代入し、この汎関数を停留させる𝑢𝑢𝑖𝑖 を求める。空間の離散化に有限要素(𝑝𝑝)法を用いると、Eq. (6)を停留させる解𝑢𝑢𝑖𝑖 の方程式は、

𝑛𝑛𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃
𝐾𝐾𝐾𝐾∑𝑛𝑛𝑝𝑝=0∑𝑝𝑝=0〈Ψ𝑝𝑝 Ψ𝑞𝑞 〉𝐾𝐾 (0) 𝑢𝑢(𝑝𝑝) + ∑𝑛𝑛𝑘𝑘=1〈𝜉𝜉𝑘𝑘 Ψ𝑝𝑝 Ψ𝑞𝑞 〉𝐾𝐾 (𝑘𝑘) 𝑢𝑢(𝑝𝑝) = 〈Ψ𝑞𝑞 〉 𝑓𝑓

(𝑞𝑞 = 0,1, … , 𝑛𝑛𝑃𝑃𝑃𝑃 ),

(7)

(𝑝𝑝)

となる。ここで、⟨ ⟩は確率変数の期待値、𝐾𝐾 (𝑘𝑘) は𝐸𝐸 (𝑘𝑘) を用いて構築された剛性行列、𝑢𝑢𝑖𝑖 は離散化された

𝑢𝑢𝑖𝑖 である。

𝑛𝑛𝑃𝑃𝑃𝑃 は𝑛𝑛𝐾𝐾𝐾𝐾 及び PC 展開次数によって定まる。本研究では確率分布が十分収束することを目指し PC 展

開次数を 3 としており、また、適用例において𝐾𝐾と𝐺𝐺に曖昧さがをもつ層𝑉𝑉𝑖𝑖 が 5 層あるため𝑛𝑛𝐾𝐾𝐾𝐾 = 10となることから、𝑛𝑛𝑃𝑃𝑃𝑃 = 285となる。モデルの節点数が𝑛𝑛で節点あたりの自由度が 3 であるため、Eq. (7)の自。このように計算
由度は3𝑛𝑛(𝑛𝑛𝑃𝑃𝑃𝑃 + 1)という大規模問題となる（e.g., 適用例では、32.0 兆自由度となる）

規模が大きくなり過ぎること、及び、Eq. (7)における左辺項の計算が煩雑となることが、このような超大規模確率有限要素解析の実現を阻んできた。本研究では、Eq. (7)における左辺項部分に関するコードを自動生成する auto code generator を開発するとともに、以下に示す主要カーネル周り及び通信に関する開発により超大規模確率有限要素解析を実現した。なお、小規模問題におけるモンテカルロシミュレーシ
ョンによる結果と本手法の結果を図 2.4.に示す。モンテカルロシミュレーションの結果は 10000 ケース
程度でようやく収束しつつあり、その収束先が本手法による 1 回の解析により得られていることが分かる。
①-(c)-1). ソルバーアルゴリズム
確率有限要素解析のほぼすべてのコストは Eq. (7)の求解にかかるため、この連立方程式を求解するソルバーを開発する。超大規模問題となるため、「富岳」等を使っても直接法で計算することは現実的でなく、メモリ削減型の前処理を使った反復法ソルバーを使うこととなる。以後、説明を明瞭にするため、求解対象の方程式 Eq. (7)を確率基底方向に書き下し- 45 -

𝐴𝐴𝑢𝑢 = 𝑏𝑏

(8)と表記する。ここで、𝐴𝐴はサイズ3𝑛𝑛(𝑛𝑛𝑃𝑃𝑃𝑃 + 1) × 3𝑛𝑛(𝑛𝑛𝑃𝑃𝑃𝑃 + 1)の正定値対称疎行列、𝑢𝑢, 𝑏𝑏はサイズ3𝑛𝑛(𝑛𝑛𝑃𝑃𝑃𝑃 + 1)の未知ベクトル・既知ベクトルとなる。性質の悪い問題の求解におけるロバスト性・収束性、前処理部のメモリ使用量、超並列環境におけるロードバランス・通信特性を考え、本研究では共役勾配法に対してマ
ルチグリッド法をベースとした精度混合型前処理を組み込んだ手法を用いる。なお、確率方向の基底の重なりよりも空間方向の基底の重なりが疎となるため、本研究では分散メモリ型並列計算においては空間
領域分割を用い、確率方向の基底はすべて同一プロセス内で保持する。以下、Algorithm 1 に沿ってソル
バーアルゴリズムを説明する。
本研究では、性質の悪い方程式もロバストに求解できる方法として、共役勾配法ベースの手法、具体的には、前処理方程式𝐴𝐴𝑧𝑧 = 𝑟𝑟の求解に他の共役勾配法ソルバーを使う、adaptive conjugate gradient methodを用いる。元の共役勾配法に関する反復を outer loop と呼び、前処理方程式の求解に用いる共役勾配法の反復を inner loop と呼ぶ(Algorithm 1 line 7-16)。前処理方程式は粗く解けばよいため、adaptive

conjugate gradient method を使うことで、精度と計算コストにトレードオフがある計算手法をソルバー
内で組み合わせることが可能となる。本研究では、通常の有限要素法における超並列ソルバー[7]、[8] を参考に、前処理においてマルチグリッド法を使うことで低次モードを効率的に求める。ここでは、四面体
一次要素により対象方程式を粗く解き(inner coarse loop (Algorithm 1 line 12))、その結果を四面体二次要素による対象方程式の初期解として用いることで、四面体二次要素による対象方程式の求解(inner ﬁneloop (Algorithm 1 line 14)) に要する反復数を削減している。これらの inner loop conjugate gradientsolvers においてはメモリ使用量が少なく超並列時においても負荷分散が保たれるブロックヤコビ前処理を用いる(節点・確率基底毎に 3×3 の局所的な前処理を適用する)。また、前処理は近似計算でよいために
inner loop は単精度で計算し、outer loop は倍精度で計算する。結果として、本手法はロードバランスに優れ、また、省メモリでありつつ通常のブロックヤコビ前処理と比べて収束性・演算量・データ移動量ともに大幅に削減されると期待される。このアルゴリズムの計算コストの大部分は inner loop において各反復で計算される単精度での行列ベクトル積となる。
そこで次にこの単精度行列ベクトル積の計算方法について説明する。

- 46 -

①-(c)-2). 計算ノード内での疎行列ベクトル積
Eq. (8)における𝐴𝐴のサイズは巨大となるものの、空間離散化に関する部分と PC 展開による確率空間の離散化部分に分離することができ、
𝑛𝑛

𝑛𝑛

𝐾𝐾𝐾𝐾
𝑃𝑃𝑃𝑃
∑𝑝𝑝=0
𝐴𝐴 = ∑𝑘𝑘=0
Ψ (𝑘𝑘,𝑝𝑝) 𝐾𝐾 (𝑘𝑘) 𝑃𝑃(𝑝𝑝)(9)と記載できる。ここで、Ψ (𝑘𝑘,𝑝𝑝) , 𝐾𝐾 (𝑘𝑘) , 𝑃𝑃(𝑝𝑝) はそれぞれ、3𝑛𝑛(𝑛𝑛𝑃𝑃𝑃𝑃 + 1) × 3𝑛𝑛, 3𝑛𝑛 × 3𝑛𝑛, 3𝑛𝑛 × 3𝑛𝑛(𝑛𝑛𝑃𝑃𝑃𝑃 + 1)の大- 47 -

きさを持つ行列である。𝑃𝑃(𝑝𝑝) は𝑢𝑢(𝑝𝑝) ⇐ 𝑃𝑃(𝑝𝑝) 𝑢𝑢として𝑢𝑢を𝑢𝑢(𝑝𝑝) にマッピングする定数行列、Ψ (𝑘𝑘,𝑝𝑝) は PC 基底に𝑇𝑇関する行列で空間方向で共通なスパースな行列となる(Ψ (𝑘𝑘,𝑝𝑝) = �𝑎𝑎0 𝐼𝐼, 𝑎𝑎1 𝐼𝐼, … 𝑎𝑎𝑛𝑛𝑃𝑃𝑃𝑃 𝐼𝐼� 、ここで𝐼𝐼は3𝑛𝑛 × 3𝑛𝑛の単位行列)。主要な計算コストは3𝑛𝑛 × 3𝑛𝑛の疎行列𝐾𝐾 (𝑘𝑘) と、𝑃𝑃(𝑝𝑝) に𝑢𝑢をかけて生成されるサイズ3𝑛𝑛 × (𝑛𝑛𝑃𝑃𝑃𝑃 + 1)の密行列𝑈𝑈 = {𝑢𝑢(0) , 𝑢𝑢(1) , … , 𝑢𝑢(𝑛𝑛𝑃𝑃𝑃𝑃 ) }の掛け算にかかることとなる。このように確率有限要素法（SFEM）における主要カーネルは共通の疎行列に多数の右辺ベクトルをかけ合わせる generalized sparse matrixvector product (GSpMV)となり、節点毎に各節点に接続する𝑚𝑚個の節点による寄与を3 × 3𝑚𝑚と3𝑚𝑚 ×(𝑛𝑛𝑃𝑃𝑃𝑃 + 1)の密行列・密行列積として計算可能となるため、通常の有限要素法の SpMV と比べて高い性能が期待される。その一方で、行列サイズが小さいために SGEMM をそのまま使うだけだと性能がでにく
く、また、Ψ (𝑘𝑘,𝑝𝑝) の項の計算が通常の GSpMV に加えて必要となるためにこのカーネルに適した計算アルゴリズムを設計する必要がある。本研究では、𝐾𝐾 (𝑘𝑘) を CRS 形式で保持した上で節点毎に GSpMV を計算し、その結果を確率方向の適切なベクトル成分に足しこむことでΨ (𝑘𝑘,𝑝𝑝) の効果を反映する。以下、CPU,GPU それぞれに適した計算アルゴリズムを示す。

CPU においては、図 2.5.の npc 方向に SIMD 並列化を施すと共に、u の L1, L2 キャッシュへのソフ
トウェアプリフェッチを行ない、各 CPU におけるレジスタ数、演算レイテンシ・キャッシュレイテンシ・
メモリレイテンシを考慮してループアンローリング展開数・ソフトウェアプリフェッチのリードループ
数を選択することで幅広い CPU アーキテクチャで高性能が期待できる。行列Ψ (𝑘𝑘,𝑝𝑝) は各行・各列ともに数個の成分しか持たない非常に疎な行列となり密行列として計算するとコストが高くなるため、計算した𝐵𝐵(𝑘𝑘) ⇐ 𝐾𝐾 (𝑘𝑘) 𝑈𝑈 の値(図中 blocal)non-SIMD 計算で操作した上で global な結果配列𝑏𝑏に足しこむ。

GPU においては、図 2.5.の節点方向𝑛𝑛のループと npc 方向のループを collapse したループをスレッド並列化することでループ長を増加させレイテンシを隠蔽するとともに、スレッド内でのデータアクセスの coalesed access 化によりメモリアクセス効率を改善している。GPU においては、レジスタ数の制限
から CPU 同様にテンポラリ配列を用いて𝐵𝐵(𝑘𝑘) の値を操作しようとするとレジスタスピル・フィルが多発

し性能が著しく低下するため、Ψ (𝑘𝑘,𝑝𝑝) の効果を反映した結果をグローバルな結果配列𝑏𝑏に対して atomic

add で直接足しこむ。

図 2.5. 行列-ベクトル積の擬似コード
実際の計算では図中の 1 次元問題の代わりに 3 次元で実施した。

- 48 -

①-(c)-3). 計算ノード間マトリックス・ベクトル積の通信アルゴリズム
行列ベクトル積にかかる計算量は領域の体積に比例する一方で、通信量は領域の表面積に比例する。そのため、メモリ容量当たりの有限要素法の節点数𝑛𝑛が通常の deterministic な有限要素法と比べて小さくなる SFEM においては、通信がボトルネックになりやすい。そこで本研究では、行列ベクトル積における袖通信コストの削減アルゴリズム、及び、ネットワーク特性に適した問題分割・プロセス配置方法を開発する。
まず、前処理において使われる FP32 の通信変数の代わりに 16 bit integer を使うことで通信量を削減する。ここでは、 ﬂoat 変数 32 個を 1 セットとし、ﬂoat ﬂoatarray[32]の代わりに unsigned shortintarray[32], ﬂoat maxval = max(ﬂoatarray[:]), minval = min(ﬂoatarray[:]) を送受信に用いる。すなわ
ち、in-tarray[:] = (ﬂoatarray[:]-minval)/(maxval-minval)*65,535 として送信し、受信側で ﬂoatarray[:]= (intarray/65535)*(maxval-minval) + minval で展開する。これにより、メッセージサイズが(2 × 32 +4 × 2)/(4 × 32) = 1/1.777 に削減される。隣接通信時において、同じ PC 成分において隣接する節点値のオーダーは同程度になると予想されるため、隣接する節点 32 個をまとめて送受信することで精度劣化を抑える。同時に、MPI の persistent communication を利用することで、通信・計算のオーバーラップに対応したシステムにおいて通信・計算のオーバーラップを促進する。MPI 領域間で共有される節点に対する行列ベクトル積を先に計算しておき、その後、共有節点の通信と、MPI 領域内部の節点の行列ベク
トル積を並行して実行する。
大規模システムにおいては計算ノード間の通信バンド幅・レイテンシは均質ではないことが一般的であるため、ネットワークトポロジを意識した問題分割・プロセス配置が通信性能向上のために重要となる。例えば Sunway TaihuLight において採用されている four-stage tapered fat-tree network, Piz Daint
において採用されている Dragonﬂy-based network では、injection bandwidth と bisection bandwidthにそれぞれ約 7 倍と 2 倍の性能差があり、また本研究において計算に用いる「富岳」における 6D torus
network においてはトーラス内のプロセス配置により通信ホップ数・通信パスの輻輳度が変化することとなる。そこで、本研究ではネットワークが 3 階層のレイテンシ・バンド幅特性からなると仮定した、空間領域分割・プロセス配置方法を開発した(図 2.6.)。まずは大域的通信の削減のため、有限要素モデルの薄い方向(この場合は x−y 方向) を 2D sliced grid で分割する。3 次元トーラスにおいては X−Y 軸に沿ってプロセスを配置し、tree 型ネットワークでは各スイッチが一つのブロックに対応するようにプロセスを配置する。続いて、グラフ分割ソフトウェアを使って各ブロックをトーラスの Z 方向ないしスイッチ
内の計算ノード数に分割し、3 次元トーラスにおいてはこれを z 方向に並び替えてトーラスの Z 軸に沿ってプロセスを配置する。最後に、グラフ分割ソフトウェアを使って計算ノード内の複数プロセス用に領域を分割する。これにより、3 次元トーラスネットワークにおいては通信ホップ数の削減と通信パスの競合を抑えることが可能となり、tree 型ネットワークにおいては経由するスイッチ数と大域的な通信量の削減が可能となる。
①-(d). 性能の計測方法①-(d)-1). 性能計測で用いたアプリケーション
図 2.2.の南海領域の地殻モデルを模擬した 11 層の地殻モデルを使って性能計測を実施する。問題特性が南海モデルと類似したものとなるよう、南海モデルと同一の物性・要素サイズ・メッシュ生成要件を用- 49 -

いる(図 2.3.)。結果、層番号 1～5 の𝐾𝐾, 𝐺𝐺 が確率変数となり、層番号 6～11 は deterministic な物性を持つこととなる。この領域の中心に点震源を与えた際の変位の確率応答を相対誤差𝜖𝜖 = 10−8で求める。生成した有限要素モデル及び震源を𝑥𝑥, 𝑦𝑦方向に複製して配置することで計算ノード当たりの問題規模が一定の weak scaling 問題セットを作成する(表 2.1.)。静的な問題ではグリーン関数の影響域が広く問題規模の拡大とともに反復法ソルバーの収束特性が悪くなるため、共役勾配法一反復当たりの実行時間とFLOPS のピーク性能比を用いて weak scaling モデル間での性能を比較する。

図 2.6. メッシュネットワークトポロジーに対する階層的な領域分割・プロセス配置方法
ノード形状(NX, NY, NZ)=(12, 12, 4)で 1 ノードあたり 4 プロセスの場合の例

上記問題において、提案手法の性能を PCGE SFEM ソルバー、及び、deterministic な有限要素法と比較する。PCGE SFEM ソルバーは、3 × 3 ブロックヤコビ前処理を用いた共役勾配法ソルバーであり、提案手法から精度可変型のマルチグリッド前処理を省いたものに相当する。Deterministic な有限要素法
としては、Tensor Core などの特殊な計算機構を前提としない地殻変動解析の SoA である[24], [23]を「富岳」用にチューニングしたものと比較する(以後 FEM ソルバーと表記する)。提案手法および FEM ソル
バーにおいては、前処理部は単精度で、それ以外は倍精度で計算される一方、PCGE SFEM ソルバーは
すべて倍精度で計算される。PCGE SFEM・FEM ソルバー共に行列ベクトル積における integer 通信と通信オーバーラップは無効化してある。計算時間のうち、ソルバー部にかかった時間を計測対象とする。
なお、IO コストの割合が大きくなると考えられる全系実行時の適用例においても、ソルバーの実行時間と IO 時間の比は 1885.4 s:115.9 s = 16.2:1 に抑えられており、ソルバー部のみの計測でアプリケーショ
ン全体の特性をとらえることができる。
①-(d)-2). 性能測定環境
開発手法は幅広いシステムで性能が得られると期待されるが研究では最大規模の総メモリ容量を搭載し南海モデルの計算が実施可能となる「富岳」にて性能計測を実施する。また、近年幅広く GPU シス- 50 -

テムが活用されていることを念頭に、NVIDIA GPU におけるカーネルレベルの性能を測ることで提案手法の performance portability を確認する。なお、GPU 環境としては、NVIDIA A100 PCIe 80 GB GPUを使用する。A100 PCIe 80 GB GPU は FP64、FP32 の理論性能がそれぞれ 9.7 TFLOPS、19.5 TFLOPSとなっており、GPU メモリ帯域幅は 1935 GB/s となっている。

表 2.1．性能計測・適用例モデルの構成
モデルセット W と S は、weak scaling モデルセットと strong scaling モデルセットを示す。全モデルにおいて𝑛𝑛𝑃𝑃𝑃𝑃 = 285, 𝑛𝑛𝐾𝐾𝐾𝐾 = 10を用いた。

①-(e). 性能計測結果①-(e)-1). カーネル性能
まず主要計算である単精度行列ベクトル積カーネルの性能を測った。
「富岳」の A64FX CPU 上では、
3.07 TFLOPS FP64 peak の 48.54%に相当する 1.49 TFLOPS, 1024 GB/s のピークメモリバンド幅の24.14%に相当する 247 GB/s が得られた一方、A100 PCIe 80 GB 上では、9.7 TFLOPS FP64 peak の44.8%に相当する 4.37 TFLOPS, 1935 GB/s のピークメモリバンド幅の 29.1%に相当する 562 GB/s が得られた。このように、通常の有限要素法で使われる SpMV と比べて飛躍的に高い性能が CPU/GPU の双方で得られている。
①-(e)-2). Time-to-solution
提案手法の time-to-solution を表 2.1.記載の model W-2 (「富岳」4608 node)にて計測した(図 2.7.)。
まず参考値として、「富岳」用にチューニングされた deterministic な FEM ソルバーによりモンテカル
ロシミュレーションを実施した場合の実行時間を見積もる。提案手法と同じ計算ノード数を用いて FEM
1 ケース計算した際の実行時間は 6.07 s であったため、確率応答分布が収束すると期待される 10000 ケースの計算には 60707 s かかることとなる。これに対して、開発した超並列計算用の SFEM アルゴリズ- 51 -

ムにより、前処理に工夫をしない状態(PCGE) でも 2051 s で同等の計算が可能となっている。さらに精度混合演算を活用したマルチグリッド前処理により、通信回りに工夫をしない状態でも計算時間は 216.7s に改善している(FP32 通信、通信オーバーラップ無し)。通信を integer 化しても反復数はほとんど変わらなかったため、通信量がメッセージ長の減少分だけ減り、実行時間は 189.3 s に短縮された。また、
通信オーバーラップにより計算時間はさらに短縮され、最終的には deterministic 解析を使った場合の332 倍速、PCGE の 11.2 倍速に相当する 182.6 s に短縮された。ここでの FP64 ピーク性能比は 20.4%、
メモリバンド幅の活用率は 28.1%と非構造格子有限要素法としては非常に高い性能を実現している。

N/A

N/A

図 2.7. 「富岳」上でのモデル W-2、W-5、W-7 におけるソルバーの性能比較

ノード数の増加に伴う提案手法の有効性を確認するため、ノード当たり自由度を一定にしてよりノード数の多い W-5 (36864 node) と W-7 (147456 node) で同様の計測を実施した(図 2.7.)。なお、モデルサ
イズの増大に伴い連立方程式の収束性も劣化し反復数が増加することとなるため、実行時間を比較するのでなくノード当たり FLOPS 値でモデル間の性能を比較する。Integer 通信による通信量の削減、および、通信オーバーラップの効果は全モデルで効果がでており、特にノード数が増加するとその効果が顕著- 52 -

になっていることがわかる。また、36864 node においては適切なプロセス配置をすることが高い性能を得るために重要であることがわかる。結果、36864 node 時でピーク性能比 20.8%、147456 node 時でピーク性能比 19.8%と、4608 node でのピーク性能比 20.4%とほぼ同等の性能を維持できていることがわかる。

図 2.8. 「富岳」における weak scaling①-(e)-3). スケーラビリティ
次に、提案手法の「富岳」全系までの scalability を計測した(図 2.8.)。問題規模の増加に伴い実行時間は増えているものの、対象問題の収束特性が劣化しているためであり、ソルバーのピーク性能比で見ると2304 node では FP64 peak 22.6%となっていたのに対し、全系 147456 node では FP64 peak 19.8%と高い性能が維持されており、FLOPS 値でみた weak scaling efﬁciency が 87.7%となっている。全系実行時のピーク性能は 89.8 PFLOPS であり、32 兆自由度の巨大な線形連立方程式を 233.8 s という短時間で計算できている。なお、ネットワークトポロジを考慮したプロセス配置をしない場合においては、全系実行時に性能が大幅に劣化していたことから、「富岳」全系で高い性能を得るためにはネットワークト
ポロジを意識したプロセス配置が必須であることがわかる。
①-(e)-4). 適用例における性能
最後に南海領域の適用例計算を「富岳」全系で実行する際のアプリケーション性能を計測した。通常、モデル分割では METIS などのグラフ分割手法によりモデルを compute node 数で分割した上で、各領域を node 内のプロセス数に分割する(以後、2-stage METIS と呼ぶ)ことが多い。この方法は表 2.2.に示
すように、要素数や節点数のインバランスを抑えつつ通信節点数を抑えることができるものの、3 次元トーラスにおいて適切にプロセスを配置できないため、対象とする node-shape (Nx, Ny, Nz) = (48, 8, 48)にプロセスを配置した際の隣接通信時の通信ホップ数の平均値が 21.5 hops、最大ホップ数が 79 hops と非常に大きくなる。一方で、提案する partitioning 手法を使った場合、2-stage METIS 同様の節点数・要素数インバランスや通信量を維持しつつ、隣接通信時の通信ホップ数の平均値が 3.57 hops、最大値が12 hops に減り、通信性能改善が期待される。表 2.2.に計測した性能を示す。先に述べた性能計測問題での性能と比べ、2-stage METIS を使った場合では FP64 ピーク性能比が著しく低下した一方で、提案手- 53 -

法を用いることで性能劣化が抑えられていることがわかる(FP64 ピーク性能比 19.1%)。これにより、複雑なメッシュ構造を持つ 32 兆自由度の大規模問題を 878 s という短時間で求解可能となっている。

表 2.2. 適用例問題に対する分割の性質、及び、性能計測結果
分割の性質は MPI プロセスあたりの値を示している。

①-(f). まとめと今後の展望
本研究は、「富岳」全系までスケールし、かつ、高速な、Uncertainty Quantiﬁcation を可能とする超大規模自由度確率有限要素解析を開発した。ここでは非構造四面体二次要素で確率空間・実空間を離散化し、精度混合演算 implicit iterative solver により倍精度で求解している。SoA と目される手法と比較して 224 倍の高速化を達成するとともに、実問題の適用例を解くことでその有効性を示した。地震分野の被害軽減や現象理解への貢献が期待されるとともに、同様に Uncertainty Quantiﬁcation が必要とされる他分野への波及効果も期待される。
また、本研究は大規模問題における Uncertainty Quantiﬁcation において、「富岳」のような大規模ス
パコン上での超大規模な capability computing が、小規模スパコンを多数使った capacity computing を凌駕しうる可能性を示した。従来の大規模スパコンでは本研究のような解析を行うには少し小さかった
が、exascale を見据える規模になって初めて実現可能となった。より大規模なスパコンが開発されることにより、このような大規模スパコン上における大規模 capability computing が一般的となり、Uncertainty Quantiﬁcation やこれに類する解析がより高速化・高詳細化されていくことが期待される。
本研究で示した手法の開発により、従来から考慮されてきた観測誤差の影響に加えて、モデル誤差である３次元的に複雑な地下構造の物性の曖昧さまでをも含めた、大地震発生時やその後の afterslip、あるいは自発的に発生する slow slip の推定によるプレート境界断層すべりの現状把握における Uncertainty
Quantiﬁcation の定量的分析が可能となると期待される。定量的な Uncertainty Quantiﬁcation を伴った現状把握にもとづいて、断層すべりによる応力変化が周囲の断層にもたらす影響、さらにはその後のすべり発展の推移予測もできるようになると期待され、将来の「地震予報システム」の必須技術となると期待される。

- 54 -

[1] United States Geological Survey, Earthquakes(https://www.usgs.gov/programs/earthquake-hazards/earthquakes).
[2] T.H. Jordan, SC15 Invited Talk Spotlight: Societal Impact of Earthquake Simulations at Extreme
Scale.
[3] 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
SPECFEM3DGLOBE, Proceedings of the International Conference on High Performance Computing,
Networking, Storage and Analysis (SC’08). IEEE Press, 2008, pp. 1–11.
[4] 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.” Proceed-ings of theInternational Conference on High Performance Computing, Networking, Storage and Analysis(SC’14). IEEE Press, 2014, pp. 3–14.
[5] 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-Pﬂops 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.
[6] 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 the International Conference for High Performance Computing, Networking,Storage and Analysis, pp. 1–16, 2017.
[7] T. Ichimura, K. Fujita, S. Tanaka, M. Hori, M. Lalith, Y. Shizawa, and H. Kobayashi. “Physicsbased urban earthquake simulation enhanced by 10.7 BlnDOF x 30 K time-step unstructured FE nonlinear seismic wave simulation.” Proceedings of the International Conference on High Performance
Computing, Networking, Storage and Analysis (SC’14). IEEE Press, 2014, pp. 15–26.
[8] 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 to Enhance Comprehensive Earthquake Simulation.” Proceedings ofthe International Conference on High Performance Computing, Networking, Storage and Analysis(SC’15). ACM, 2015, p. 4.
[9] 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 ﬁnite elements with artiﬁcial intelligence and transprecision computing.” Proceedingsof the International Conference on High Performance Computing, Networking, Storage and Analysis(SC’18). ACM, 2018, Article No. 49,[10] Matthew E. Pritchard, Richard M. Allen, Thorsten W. Becker, Mark D.Behn, Emily E. Brodsky,
Roland Burgmann, Cindy Ebinger, Jeff T. Freymueller, Matt Gerstenberger, Bruce Haines, Yoshihiro- 55 -

Kaneko, Steve D. Jacobsen, Nate Lindsey, Jeff J. McGuire, Morgan Page, Sergio Ruiz, Maya Tolstoy,
Laura Wallace, William R. Walter, William Wilcock, and Harold Vincent, New Opportunities to Study
Earthquake Precursors, Seismol. Res. Lett. XX, 1?4, doi: 10.1785/0220200089, 2020.
[11] H. Noda and T. Hori, Under what circumstances does a seismogenic patch produce aseismictransients

in

the

later

interseismic

period?,Geophys.

Res.

Lett.,

41,

7477-7484,

doi:10.1002/2014GL061676, 2014.
[12] Y. Kaneko, L. M. Wallace, I. J. Hamling, and M. C. Gerstenberger, Simple physical model for the
probability of a subduction-zone earthquake following slow slip events and earthquakes: Application
to the Hikurangi megathrust, New Zealand, Geophys. Res. Lett. 45, 3932-3941, doi:
10.1029/2018GL077641, 2018.
[13] M. C. Gerstenberger, Y. Kaneko, B. Fry, L. Wallace, D. Rhoades, A. Christophersen, and C.
Williams, Probabilities of large earthquakes in central New Zealand (p. 24). Wellington: Published
Report to New Zealand Natural Hazards Platform, 2017.
[14] C. A Williams and L.M. Wallace, The impact of realistic elastic properties on inversions of shallow
subduction interface slow slip events using seaﬂoor geodetic data, Geophys. Res. Lett., 45, 7462-7470.
https://doi.org/10.1029/2018GL078042, 2018.
[15] K. Obara and A. Kato, Connecting slow earthquakes to huge earth- quakes, Science 353 (6296),
253-257, doi: 10.1126/science.aaf1512, 2016.
[16] T. Hori, Mechanisms of separation of rupture area and variation in time interval and size of greatearthquakes along the Nankai Trough, southwest Japan. J Earth Simul 5:8-19, 2006.
[17] Takane Hori, Mamoru Hyodo, Shin’ichi Miyazaki and Yoshiyuki Kaneda, Numerical forecastingof the time interval between successive M8 earthquakes along the Nankai Trough, southwest Japan,
using ocean bottom cable network data, Mar Geophys Res, doi: 10.1007/s11001-014-9226-8, 2014.
[18] A. Nakanishi, N. Takahashi, Y. Yamamoto, T. Takahashi, S. O. Citak, T. Nakamura, K. Obana, S.
Kodaira, and Y. Kaneda, Three-dimensional plate geometry and P-wave velocity models of the
subduction zone in SW Japan: Implications for seismogenesis, Geol. Tectonics Subduction Zones
Tribut. Gaku Kimura Geol. Soc. Am. Spec. Pap., 534, 69-86, 2018.
[19] K. Koketsu, H. Miyake, Afnimar, and Y. Tanaka, A proposal for a standard procedure of modeling
3-D velocity structures and its application to the Tokyo metropolitan area, Japan. Tectonophysics, 472,
290-300, 2009.
[20] K. Koketsu, H. Miyake, and H. Suzuki, Japan integrated velocity structure model version 1, In
Proceedings of the 15th world conference on Earthquake Engineering, Lisbon, Portugal, 24–28
September, Number 1773, 2012.
[21] S. Aoi, Y. Asano, T. Kunugi, T. Kimura, K. Uehira, N. Takahashi, H. Ueda, K. Shiomi, T.
Matsumoto and H. Fujiwara, MOWLAS: NIED observation network for earthquake, tsunami andvolcano, Earth, Planets and Space, 72:126 https://doi.org/10.1186/s40623-020-01250-x, 2020.
[22] Tsuyoshi Ichimura, Ryoichiro Agata, Takane Hori, Kazuro Hira-hara, Chihiro Hashimoto, Muneo
Hori, Yukitoshi Fukahata, An Elastic/Viscoelastic Finite Element Analysis Method for Crustal- 56 -

Deformation using a 3D Island-scale High-ﬁdelity Model, Geophysical Journal International, 114-129,
206, 2016.
[23] Takuma Yamaguchi, Kohei Fujita, Tsuyoshi Ichimura, Akira Naruse, Jack C. Wells, Christopher
J. Zimmer, Tjerk P. Straatsma, Muneo Hori, Lalith Maddegedara, and Naonori Ueda. 2020. LowOrder Finite Element Solver with Small Matrix-Matrix Multiplication Accelerated by AI-Speciﬁc
Hardware for Crustal Deformation Computation. In Proceedings of the Platform for Advanced
Scientiﬁc Computing Conference (PASC ’20). Association for Computing Machinery, New York, NY,
USA, Article 16, 1–11. DOI:https://doi.org/10.1145/3394277.3401860[24] Tsuyoshi Ichimura, Kohei Fujita, Takuma Yamaguchi, Akira Naruse, Jack C. Wells, Christopher
J. Zimmer, Tjerk P. Straatsma, Takane Hori, Simone Puel, Thorsten W. Becker, Muneo Hori, andNaonori Ueda, 416-PFLOPS fast scalable implicit solver on low-ordered unstructured ﬁnite elements
accelerated by 1.10-ExaFLOPS kernel with reformulated AI-like algorithm: For equation-based
earthquake modeling, SC19 Research Poster.
[25] Ghanem, R.G. and Spanos, P.D., Stochastic ﬁnite elements, a spectral approach, Springer, Berlin,
1991.
② シミュレーション内で生成されるデータの学習による高詳細粘弾性地殻変動解析の高速化②-(a). はじめに
時刻歴の偏微分方程式(PDE)における陰解法ベースのシミュレーションにおいては、大規模疎行列の連立方程式の求解が非常に高コストとなるため、マルチグリッド法など並列計算機をスケーラブルに活用できる equation-based な求解手法が開発されてきた。その一方で、近年ではディープラーニングなどのデータ駆動型の手法も開発されており、連立方程式ソルバーの分野においても equation-based な手法とデータ駆動型の手法を組み合わせることでさらなる効率化が期待されている。大規模課題を解く場合、対象となる計算機アーキテクチャでアルゴリズムが低コストで実行できることに加え、多くの計算ノードへのスケーラビリティが必須であり、これらの要件を満たしつつデータ駆動型の手法を equation-based
ソルバーに統合する必要がある。PDE に基づく連立方程式の解法にデータ駆動型手法を用いた過去の研究では、小規模な問題に対しては[1][2]、大規模な問題に対しては [3-5]などの手法が開発されている。例
えば、[5]では、時間発展問題の PDE の一つである波動方程式を解く際に、過去の時間ステップの解を学習して次のステップの解を高精度に予測し、反復ソルバーの反復数を減らす方法を提案している。ここでは、通信不要の局所学習により局所モードを予測し、equation-based な手法であるマルチグリッド法でグローバルな誤差低減を保証することで、計算精度を落とすことなく「富岳」で効率的な並列計算を実現した。このような equation-based な手法とデータ駆動型手法を組み合わせたアルゴリズムは、幅広いア
プリケーションに応用可能と期待されるが、対象となる PDE や離散化の特性によってその有効性が変化すると考えられるため、データ駆動型の初期解推定と equation-based なソルバーの組み合わせに関する
さらなる研究が必要である。そこで本研究では、粘弾性地殻変動解析を対象に、データ駆動型の初期解予測を用いたスケーラブルな陰解法非構造化有限要素ソルバーを開発し、「富岳」の超並列計算環境での性能を測定する。このような粘弾性応答解析では、応力の収束を確保し、地震発生過程の解析や次の地震の可能性を評価するために、101－2 m 程度の非常に高解像度なメッシュで断層面付近を離散化しつつ、106m ×106－7 m×105－6 m の広い領域を対象とする必要がある。さらに、地震後の粘弾性緩和による応力－7- 57 -

状態の変化を追跡するためには、102―3 時間ステップの時刻歴解析を行う必要があり、109 自由度を持つ
非常に大きなモデルによる線形方程式系を 102―3 回解くこととなる。粘弾性解析は、これまでの研究でデータ駆動型手法が適用されてきた波動伝播問題と比較して、グリーン関数の影響範囲が広く、低次モードが支配的な問題となる。本研究では、データ駆動型手法と適切なマルチグリッドソルバーを組み合わせることで、ソルバー全体を効率化することを目指した。大地震が発生するプレート沈み込み帯の一つである
南海沖の地震後の粘弾性地殻変動の適用例で、提案手法の有効性を確認した。
②-(b). 対象問題
本研究では地殻を Maxwell 粘弾性体としてモデル化し、非構造格子有限要素法により離散化すること
で地殻の時間発展問題を求解する。ここでは、複雑な形状や不均質な材料特性を持つ領域の離散化に適す
る四面体二次要素を用いる。粘弾性地殻変動解析における時間発展については、［6］に基づき計算する。
この計算における主要コストは各時間ステップ𝑖𝑖における正定値対称な疎行列𝐾𝐾 𝑣𝑣 に関する
𝐾𝐾 𝑣𝑣 𝛿𝛿𝑢𝑢𝑖𝑖 = 𝑓𝑓 𝑖𝑖(1)の求解に係るため、このコストを削減することを考える。ここで、𝛿𝛿𝑢𝑢𝑖𝑖 は変位増分、𝑓𝑓 𝑖𝑖 は外力である。
②-(c). 開発手法
過去の時間ステップの結果を利用するデータ駆動型手法により、式(1)の効率的なソルバーを開発する。
ここでは、超並列計算機システムにおいて、均一な負荷分散を確保しつつ通信量を削減可能なソルバーを設計している。具体的には、データ駆動型手法に基づいて高精度な解を推定し、マルチグリッドベースの反復法ソルバーの初期解として用いることで、反復型マルチグリッドソルバーの反復回数の削減、及び、プログラム全体のコスト削減を目指す。
本研究では、陰解法ソルバーの反復回数削減に直結する、より良い初期解を推定するために［5］に基づくデータ駆動型予測手法を使用する。ここでは、動的モード分解(DMD)の考え方［7］を拡張し、𝑖𝑖 – 1

番目の時間ステップまでの計算結果を学習し、𝑖𝑖番目の時間ステップにおける初期解を推定する。ここでは、解析領域を小領域に分割し、各領域において、(2)に記載したようにトレンドを含む解は解析コストの低い equation-based な方法である 2 次の Adams-Bashforth 法を用いて推定し、残りのモードはデータ駆動型手法により超並列計算に適した形で局所的な時間・空間で推定する。
𝑖𝑖
𝑖𝑖
𝑖𝑖−1
= 𝛿𝛿𝑢𝑢𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎+ 𝐶𝐶(𝛿𝛿𝑢𝑢𝑖𝑖−1 − 𝛿𝛿𝑢𝑢𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎)
𝛿𝛿𝑢𝑢𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝(2)

なお、METIS［8］を用いて各 MPI プロセスに割り当てられた領域を重複しないように分割し、各小領域に含まれる節点の変位増分を同一の小領域内の時刻歴データを用いて時間発展オペレータ𝐶𝐶を推定し、初期解の推定に用いている。以上のように、小領域ごとに推定を行うため、挙動予測に必要なモード数（学習に必要なデータ数）を減らすことができ、大域計算を避けることで高効率な超並列計算が可能となる。さらに各小領域における推定計算はスレッド並列で計算できマルチコアの SIMD 計算に適した形となっている。

- 58 -

本研究では、この推定手法を前処理付き共役勾配法ソルバーの初期解として用いることで解析の高速化を図る。その際、対象となる四面体二次要素要素モデルの求解の前処理にて 3 段階のマルチグリッド
ソルバーを使用する。最終的な解の精度を保証するために外側ループは FP64 で計算し、前処理に使用するマルチグリッドソルバーは FP32 で計算することで、データ転送サイズ、通信サイズ、メモリ使用量を半分に削減した。
②-(d). 性能測定条件
実際の地殻変動問題を模擬して作成した性能測定問題を用いて性能を計測した(表 2.3.)。開発した手法の有効性を示すため、2 次の Adams-Bashforth 法にて初期解を設定した 3×3 ブロックヤコビ前処理による共役勾配法ソルバー（PCGE）、及び、マルチグリッド前処理による共役勾配ソルバー（以下マルチ
グリッドソルバーと呼ぶ）と性能を比較した。また、開発した初期解の予測手法の効果を測定するため、Adams-Bashforth 法により初期解を設定した PCGE と、開発したデータ駆動型手法により初期解を設定した PCGE の性能を比較した。初期段階での粘弾性応答の揺らぎが抑制され、粘弾性解析全体のほとん
どの時間ステップの特性に類似している時間ステップ 21-30 で、経過時間と反復回数を測定した。計算には MPI と OpenMP によるハイブリッド並列プログラミングモデルを使用した。理化学研究所計算科
学研究センターのスーパーコンピュータ「富岳」で性能測定を行った。

表 2.3. 性能測定に使用するモデル
モデルセット W は weak scaling モデルセットで、モデルセット A は適用例モデルセットである。2 次の Adams-Bashforth 法により初期解を設定した PCGE ソルバーにより各モデルを求解した際の反復数も記載している。

②-(e). 性能測定結果
まず、モデル W-1 における PCGE ソルバーの反復数削減効果を確認した（図 2.9.a））。一般的に使われる二次の Adams Bashforth 法（タイムステップあたり 16377 反復）を初期解に使うことで、前ステップの解を初期解に使う場合（22807 反復）よりも反復数が削減されている。その一方で、Adams Bashforth
法の次数を 3 次に増加した場合はかえって反復数が増加しており（約 18000 反復）、ポイントワイズの予測において推定次数を増やすことが必ずしも初期解の特性改善につながるわけではないことがわかる。
また、初期解としてランダムな値を使用した場合、反復数は 15266 回とやや減る結果となった。提案するデータ駆動型手法(s=16)を PCGE の初期解として用いた場合、反復数が 5062 回に減少した。表 2.4.に示すように、ソルバー全体と比較してデータ駆動型手法による初期解推定にかかる時間はごくわずかで- 59 -

あり、二次の Adams-Bashforth 法を用いた PCGE ソルバーと比較して 3.19 倍の高速化を達成した。

図 2.9. モデル W-1（576 ノード）における各ソルバー・初期解推定手法における収束履歴
各線は、粘弾性解析の 21-30 タイムステップにおける各タイムステップの誤差収束履歴を示す。

表 2.4. Weak scaling モデル W-1 (576 ノード) における各ソルバーの経過時間とタイムステップの収束に要した反復数- 60 -

次に、モデル W-1 におけるマルチグリッドソルバーの反復回数削減の効果をみる（表 2.4.）。マルチグ
リッドソルバーにおいては粗いグリッドにより低次モードの誤差を少ない反復数で減らすことができる
ため、Adams-Bashforth 法を初期解推定に用いた場合でも総反復回数は 43+90+956+7861=8950 回とPCGE ソルバーの場合よりも少なくなる。この高効率なマルチグリッドソルバーの初期解にデータ駆動型手法を用いることで、反復数は 26 + 59 + 175 + 1893 = 2153 回に減少した（図 2.9.b））。これにより、
Adams-Bashforth 法による初期解を用いたマルチグリッドソルバーと比較して 2.68 倍の高速化を実現した。
表 2.5.に weak scaling の測定結果をまとめる。データ駆動型手法による初期解推定にかかる時間はノード数によらず一定、かつ、多数ノードにおいても良好な反復数削減効果が維持されていることがわかる。576 ノードから 73728 ノードまでの weak scaling 効率は 78.5%であり、通信削減・負荷分散特性に優れた初期解予測手法およびマルチグリッドソルバーが有効に働いていることがわかる。結果、73728 ノードの最大問題において 6.88%の高いピーク性能を達成し、Adams-Bashforth 法により初期解を設定したマルチグリッドソルバー比で 5.12 倍、PCGE ソルバー比で 196 倍の高速化を達成した。

表 2.5. Weak scaling モデルセットにおける経過時間およびタイムステップごとに必要な反復数
PCGE（W-7、W-8、*印）は 21 番目のタイムステップの値を、他のケースは 21-30 番目のタイムステッ
プの平均値を示す。

②-(f). 適用例
詳細な三次元地殻構造モデル上においてもこの提案手法が有効であることを確認するため、大地震が発生するプレート沈み込み帯の一つである日本の南海沖での地震発生後を想定し、その粘弾性地殻変動の計算を行った(図 2.10.a))。断層面の応力とｆ歪を正確に計算するため、断層近傍では 250 m、断層から遠く離れた硬い材料特性の領域では最大 4000 m の要素によりメッシュを生成している(A-1)。
さらに、得られた解析結果の収束を確認するためにメッシュ全体を一段階再分割し、解像度を倍にしたケース（A2、最小要素サイズ 125 m）も用意した（図 2.10.ｂ））。これらの有限要素モデルにより大規模な地震を想定した仮想すべり［9］に対する粘弾性応答を計算した。地震後の地表の変形を図 2.10.e)に表す。海側
では大きな陸側変形が見られ、上下方向では海溝軸付近で隆起、海溝軸の奥で沈降する典型的な逆断層地- 61 -

震の変形パターンが見られる。このような現実的な問題設定の計算における提案手法のピーク性能は8.16%であり、性能測定問題で得られた性能とほぼ同じであった。また、モデル A-1 においては提案手法
により反復数がマルチグリッドソルバーと比較して 4.2 分の 1 に削減され、結果、マルチグリッドソル
バー比で 2.9 倍、PCGE 比で 84 倍の高速化を達成した。モデル A-2 においても Adams-Bashforth 法を用いたマルチグリッド法から 3.4 倍の高速化を達成した。このように開発手法は単純化された性能測定問題だけでなく、複雑な形状やメッシュ構造を持つ適用例問題においても有効であることが分かる。

図 2.10. 適用例 a) 生成された有限要素モデル、b) 有限要素モデルの拡大図、
c) d の線分 AB の垂直断面に沿った構造、d) 入力すべり分布、e) 表面における粘弾性変位応答。

②-(g). まとめ
本研究では、データ駆動型手法による高速・スケーラブルな非構造格子有限要素法ベースの粘弾性地殻
変動解析手法を開発した。開発手法により、標準的な初期解予測手法を用いたマルチグリッドソルバーと比べて求解に必要な反復数・実行時間が削減された。開発したデータ駆動型の初期解予測手法は PDE に基づく一般的な陰的時間積分による時刻歴シミュレーションに対して適用可能であり、他の問題においても有効であると期待される。

[1] J. Sirignano and K. Spiliopoulos, “Dgm: A deep learning algorithm for solving partial differential
equations,” Journal of computational physics, vol. 375, pp. 1339–1364, 2018.
[2] J. Han, A. Jentzen, and W. E, “Solving high-dimensional partial differential equations using deeplearning,” Proceedings of the National Academy of Sciences, vol. 115, no. 34, pp. 8505–8510, 2018.
[3] T. Ichimura, K. Fujita, T. Yamaguchi, A. Naruse, J. C. Wells, T. C. Schulthess, T. P. Straatsma,
C. J. Zimmer, M. Martinasso, K. Nakajima et al., “A fast scalable implicit solver for nonlinear timeevolution earthquake city problem on low-ordered unstructured finite elements with artificial
intelligence and transprecision computing,” in SC18:International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE,- 62 -

2018, pp. 627–637.
[4] T. Ichimura, K. Fujita, M. Hori, L. Maddegedara, N. Ueda, and Y. Kikuchi, “A fast scalable
iterative implicit solver with green ’s function-based neural networks,” in 2020 IEEE/ACM 11th
Workshop on Latest Advances in Scalable Algorithms for Large-Scale Systems (ScalA). IEEE, 2020,
pp. 61–68.
[5] T. Ichimura, K. Fujita, K. Koyama, R. Kusakabe, Y. Kikuchi, T. Hori, M. Hori, L. Maddegedara, N.
Ohi, T. Nishiki et al., “152k-computernode parallel scalable implicit solver for dynamic nonlinear
earthquake simulation,” in International Conference on High Performance Computing in AsiaPacific Region, 2022, pp. 18–29.
[6] T. Ichimura, R. Agata, T. Hori, K. Hirahara, C. Hashimoto, M. Hori,and Y. Fukahata, “An elastic/viscoelastic finite element analysis method for crustal deformation
using a 3-d island-scale high-fidelity model,” Geophysical Journal International, vol. 206, no. 1, pp.
114–129, 2016.
[7] J. N. Kutz, S. L. Brunton, B. W. Brunton, and J. L. Proctor, Dynamic mode decomposition: datadriven modeling of complex systems. SIAM, 2016.
[8] G. Karypis and V. Kumar, “A fast and high quality multilevel scheme for partitioning irregulargraphs,” SIAM Journal on scientific Computing, vol. 20, no. 1, pp. 359–392, 1998.
[9] M. Hyodo and T. Hori, “Re-examination of possible great interpolate earthquake scenarios in
the Nankai Trough, Southwest Japan, based on recent findings and numerical simulations,”
Tectonophysics, vol. 600, pp. 175–186, 2013.

③ 「富岳」で開発した高効率行列ベクトル積手法の汎用化のための開発③-(a). はじめに
本プロジェクトでは地震シミュレーションにおける有限要素解析の高速化を目的に、「富岳」で高速に疎行列ベクトル積を実行できる方法を開発してきた。有限要素法における疎行列ベクトル積を「富岳」だけでなく近年普及の進む GPU 計算機においても高速実行できるように開発することで、「富岳」用に開発した高速ソルバーアルゴリズムを汎用化できると考えられる。そうすれば、企業が実務や研究開発において、より幅広い計算機環境で我々のコードを活用できると期待される。そこで本研究では「富岳」で開発した高効率行列ベクトル積手法の汎用化のための開発を実施した。
有限要素法をはじめとする物理シミュレーション分野においては、疎行列ベクトル積が主要計算になることが多い。このような疎行列ベクトル積を従来のマルチコア CPU システムにおいて実施する際には CRS などで疎行列を保持し、それを読み込んで行列ベクトル積を計算するカーネルを MPI とOpenMP で並列化し、実行することが一般的である。本プロジェクトで対象とする地震・地震工学分野
のみならずモノづくり分野などでは対象問題の形状をより正確にモデル化できる非構造格子で離散化を行うことが多いが、その場合、このカーネルはメモリアクセス卓越かつランダムアクセスが多く含まれるカーネルとなる。

- 63 -

近年、GPU などのアクセラレータが普及してきており、CPU に比べてメモリバンド幅が多い GPUを使用することでこのカーネルの高速化が期待される。GPU をはじめとしたアクセラレータを使用したプログラム開発においては、CUDA などの専用の言語を使う必要がありハードルが高いとされてきた。しかし、近年においては OpenACC といった指示文ベースの並列プログラミングモデルの開発・普及が進んでおり、これを使うことで従来よりも格段に低いコストでアクセラレータの活用が期待できる状況になっている。そこで本研究では、OpenACC により疎行列ベクトル積を GPU にポーティングすることで高速化した。この際、GPU が高い演算性能を有しており、相対的にメモリバンド幅が小さい
ことをふまえ、演算量は多いがメモリアクセス量が少ない Element-by-element (EBE)法を用いること
で、さらなる高速化を図った。
性能計測の結果、1 ABCI 計算ノード上で A100 GPU を使った場合、Intel Xeon CPU を使った場合と比べて、ブロック CRS (BCRS)形式で行列を持つ場合で 22 倍、EBE 法を使用した場合では 19 倍の高速化を実現できることが示された。また、GPU における BCRS/EBE の実行時間比は EBE が 3.4 倍速になった。これらより、GPU を活用した計算の有効性が示されるとともに、特に EBE において高い性能が得られたことがわかった。
本研究の後半ではこれらの知見を、液状化を考慮した地震応答解析に適用し、GPU を使用した解析の高速化を実現した。このアプリケーションの開発にあたっては、OpenACC の機能を活用し、CPU/GPU の双方で動くパフォーマンス・ポータブルなシングルコードのアプリケーション開発を行った。開発コードは、CPU 上では OpenMP を使って実装された CPU 用の既往手法とおおむね同等の性能を達成しつつ、GPU を使うことで CPU を使った場合の 42 倍の高速化を実現した。
③-(b). 手法③-(b)-1). 対象とする疎行列
本研究では、ものづくりや地震分野など複雑形状・非均質物性を持つ領域の数値解析に広く用いられる、低次の非構造格子有限要素法における行列方程式(1)を求解対象とする。

ここで、x は未知ベクトル、b は既知ベクトルであり、A は疎行列となる。有限要素法を用いて動的な物
理シミュレーションを行う場合、方程式(1)を繰り返し解くことになる。大規模問題の場合、この方程式
は並列計算機上で解かれることになるが、それには前処理付き共役勾配法などの反復法がよく用いられる。この際、係数行列 A とあるベクトル p の積の計算を繰り返し行うことになり、求解において最も計算コストの高い計算となるため、行列ベクトル積 Ap を高速に行うことが重要となる。本研究では疎行
列 A としては、応力を低コストで精度よく計算可能なためものづくりの実務などで広く用いられる四面体 2 次要素で離散化した場合(tet10)の係数行列と、精度は少し劣るがより低コストで計算可能な四面体
1 次要素で離散化した場合(tet4)の係数行列を計測・分析の対象とする。

③-(b)-2). 疎行列ベクトル積のアルゴリズム
疎行列ベクトル積 Ap を計算する際の一般的な計算方法は、各要素にて要素係数行列 Ae を生成しそ- 64 -

の成分を全体行列 A にとして足しこみ・メモリに保存しておいた上で、疎行列ベクトル積を計算する都度メモリから全体行列 Aを読み込んで p にかけ合わせることで計算を実施する。ここで Pe はローカルな節点番号をグローバルな節点番号にマッピングする行列である。疎行列の格納形式としては各種方法が提案されているが、本研究ではメモリアクセス量が少なく一般的な方法である CRS 形式のうち、ブロックサイズ 3 の block CRS(BCRS)を分析対象とする(ブロックサイズは 1 節点あたりの自由度に対応している。対象問題は 3 次元の問題であり、各節点が x, y, z 方向の 3 自由度を持っている。)。BCRS で疎行列 A を格納した場合の行列ベクトル積のアルゴリズムを図 2.11. に示す。このアルゴリズムでは行列の非零成分 val へのメモリアクセスがメモリアクセスの大部分を占め、メモリ律速の計算になる。CPU と比較してメモリバンド幅の大きな GPU を使うことで、この計算の高速化が期待される。

図 2.11. BCRS を用いて格納した A を用いた行列ベクトル積

一般に要素係数行列 Ae の生成に際しては、少量の要素特性データを入力に比較的データ量の多い Aeが生成されることとなる。したがって、CRS のように事前に全体行列を生成・保存し、行列ベクトル積の都度読み込む方式はメモリ使用量とメモリ転送を増やすことでその都度計算量の多い Ae を生成することを避け演算量を減らす操作に相当する。その一方で、計算量はかかるものの、データ量が小さい要
素特性データを入力に行列ベクトル積の都度 Ae を生成し右辺ベクトルにかけ合わせる Element-byElement 法(EBE 法）

を用いることで疎行列ベクトル積に係るメモリ使用量とメモリ転送を大幅に減らすことができる可能性がある。Pe, PeT はそれぞれ、グローバルなベクトルから要素ごとのベクトルの抜き出しと、その逆を行- 65 -

っているので、実質的な計算は、要素ごとの行列ベクトル積 Aepe (ただし、pe = Pep)を行い、それを全体のベクトルに足しこむことである(図 2.12.)。例えば Ae が複数の小行列の積として表せる場合など、
Ae 自体を計算し pe にかけ合わせる方法よりも、要素行列ベクトル積 Aepe を直接計算する方が計算量を低減可能となることが多いため、実際の EBE の開発にあたってはこの Aepe の計算コストを減らすような定式化・実装がポイントとなる。例えば本研究の③―(c)で性能計測の対象とする問題(tet10, tet4)
における各変数のメモリ使用量とランダムアクセス量を表 2.6.に示す。メモリ使用量は行列ベクトル積の計算において少なくとも 1 回はメモリに読み書きする必要のあるデータの量でもあるため、総メモリ
アクセス量の下限値となる。一部の変数(例えば、図 2.12.の x, Ax)は複数回アクセスされる可能性があるため、ランダムアクセスの量はメモリ使用量より多くなることもある（2 回目以降のアクセスはキャッシュへのアクセスになることも多いため、メモリへのアクセス量はこの値よりも少なくなる)。Tet10(要素行列 Ae の一辺の大きさ：10 節点× 3 自由度=30 自由度) の場合、EBE を使用することでメモリ使用量は 1/17 となりランダムアクセスも 38%減る。演算量は BCRS よりも増加すると予想されるが、GPU など近年のアーキテクチャでは、メモリバンド幅と比較して演算性能が高いため、演算量の増加
は計算時間にそれほど影響を及ぼさず、メモリアクセス量・ランダムアクセス量の削減によって計算時間が短縮されることが期待される。Tet4 (要素行列 Ae の一辺の大きさ：4 節点× 3 自由度=12 自由度)の場合、EBE を使用することで BCRS と比較してメモリ使用量は 1/4 に減っているもののランダムメ
モリアクセスは約 5 倍になっている。Tet10 ほどではないものの Tet4 においては、ランダムメモリア
クセスの増加によるコスト増よりも総メモリ転送量の低減によるコスト減が優位に働く計算機においては、疎行列ベクトル積の高速化に EBE が有効である可能性がある。また、EBE では A を保持しないた
め、メモリ使用量が少ない。そのため、より大規模な計算が可能となるなど計算速度以外の観点から
EBE の方が良い場面・アプリケーションもあると考えられる。

図 2.12. EBE のアルゴリズム

表 2.6. 行列ベクトル積におけるメモリ使用量とランダムアクセス量- 66 -

③-(b)-3). 本研究で実装したカーネル
本研究で使用したカーネルを表 2.7.に示す。アルゴリズムとしては、BCRS によって格納した疎行列を使用するもの(以後 BCRS と呼ぶ) と EBE 法を使用するもの(以後、EBE と呼ぶ)の 2 つを使用した。
BCRS: BCRS(Fugaku)-OpenMP は BCRS を「富岳」用にチューニングしたカーネルである。図 2.11.
の i1 のループを OpenMP を使ってスレッド並列化してある。さらにこの際、OpenMP スレッド間のロードバランスが保たれるように事前に各 OpenMP スレッドが担当する範囲を計算している。また、SIMD
を有効利用するために配列 val が j2 方向にメモリ上に連続になるようにし、k1, k2 のループをアンローリングした。さらに、val の L1 キャッシュへのソフトウエアプリフェッチを行うとともに、セクタキャ
ッシュを使用して、x が L1 キャッシュ上に残りやすくなるようにした。BCRS-OpenMP と BCRSOpenACC は、主に GPU 上での使用を想定したカーネルであり、この 2 つはディレクティブ以外は全く同じコードである。OpenMP/OpenACC の指示文を使って i1 のループを並列化している。BCRScuSPARSE は、GPU 上での BCRS の性能を確かめるためのレファレンスとなるカーネルで、cuSPARSE
ライブラリの疎行列ベクトル積関数 cusparseSbsrmv を使用して疎行列ベクトル積の計算を行う。この関数は BCRS 形式で格納された疎行列とベクトルの積を計算する関数である。
EBE: EBE(Fugaku)-OpenMP は既往研究で開発された「富岳」用にチューニングしたカーネルである。図 2.12.の ie のループを OpenMP を使ってスレッド並列化した。この際、カラーリングを使用することで Ax(·)+=Axe(·)における OpenMP スレッド間でのデータ競合を回避している。そして、SIMD
命令を有効利用するために、ランダムメモリアクセス部と要素ごとの行列ベクトル積の計算部を分離している。また、「富岳」の限られたレジスタ数を有効活用するため、要素ごとの行列ベクトル積の計算に関するループを分割することで、性能向上を図っている。EBE(Fugaku)-OpenACC はEBE(Fugaku)-OpenMP とディレクティブ以外は全く同じコードで、OpenMP の代わりに OpenACC
を使ってスレッド並列化を行ったものである。EBE(GPU)-OpenACC は、GPU 用のカーネルである。
ie のループを OpenACC を使って GPU 上で並列計算を行うようにした。この際、カラーリングは使わず、Ax(·)+=Axe(·)の計算において atomic add を使用した。
表 2.7. 対象カーネル一覧③-(c). 性能計測③-(c)-1). 使用した計算機
AI Bridging Cloud Infrastructure (ABCI)と「富岳」を用いて性能計測を行った。
「富岳」は理化学研
究所計算科学計算センターが運用する CPU ベースのスーパーコンピュータである。
「富岳」の各計算ノ- 67 -

ードは A64FX CPU 1 台(48 コア、base frequency: 2.0 GHz)が搭載されている。ピーク性能は 3.07 FP64
TFLOPS、 6.14 FP32 TFLOPS であり、32 GiB の HBM2 (メモリバンド幅：1.03 TB/s)を搭載している。ABCI は産総研が運用する GPU ベースのスーパーコンピュータである。ABCI には 2 種類の計算ノードがあるが、そのうち計算ノード(A)を使用した。
各計算ノード(A)は 2 台の Intel Xeon Platinum 8360Y
CPU (36 コア、base frequency: 2.40 GHz)と 8 台の NVIDIA A100 GPUs for NVLink (40GiB HBM2 を搭載)から構成されている。各計算ノードの CPU の合計ピーク性能は 5.53 FP64 TFLOPS、11.06 FP32
TFLOPS であり、512 GiB の DDR4 RAM (メモリバンド幅：408 GB/s)を搭載している。各計算ノードの GPU の合計ピーク性能は 77.6 FP65 TFLOPS、156 FP32 TFLOPS であり、メモリバンド幅は 12.4
TB/s である。
③-(c)-2). 計測方法③-(b)-3). で述べた通り、BCRS, EBE について「富岳」用及び GPU 用に実装したものを計測した。計測は、「富岳」と ABCI 上で行い、ABCI 上での計測では、GPU を使う場合と CPU のみを使う
場合の計測を行った。また、CPU のマルチコアを使用したスレッド並列化には「富岳」上ではOpenMP を使用し、ABCI (CPU)上では OpenMP と OpenACC を使用した。すべての計算は FP32 で行った。疎行列 A としては、式(1)を得る際に四面体 1 次要素で離散化した場合の係数行列(Tet4)と四面体 2 次要素で離散化した場合の係数行列(Tet10)の 2 通りを計測した。Tet4, Tet10 でのモデル自由度(=vector x の長さ) は、それぞれ、7,430,679、58,574,715 である。ABCI 上で CPU のみを使った場合はこれを 1 計算ノード上で 8 MPI プロセス×9 OpenMP スレッドで計算した。ABCI 上で GPU を使う
場合は、モデルを 8 分割したものを 1 GPU 上で計算した。また、「富岳」上での計測においては、
ABCI 上での計測において 1 GPU が計算する領域を 1 計算ノード上で 4 MPI 並列で計算した。これ
は、「富岳」8 ノード (計 49.2 FP32 TFLOPS, 8.19 TB/s)、2 Xeon CPUs (計 11.06 FP32 TFLOPS,
408TB/s) 、8 A100 GPUs (計 156 FP32 TFLOPS, 12.4 TB/s)で同じ問題規模の計算をしていることに相当する。
「富岳」上での計測では、MPI_Wtime()関数により時間計測を行い、Advanced
Performance Profiler によって計測した FLOP count とメモリアクセス量をもとに FLOPS とメモリス
ループットを算出した。ABCI(GPU)上での性能計測では MPI_Wtime()関数により時間計測を行い、Nvidia Nsight Compute によって計測した FLOP count とメモリアクセス量をもとに FLOPS とメモ
リスループットを算出した。ABCI (CPU)上での計測では、MPI_Wtime()関数により時間計測を行い、ABCI(GPU)上または「富岳」上で同一カーネルを計測した時と FLOP count とメモリアクセス量が同じであると仮定し、FLOPS とメモリスループットを推定した。
③-(c)-3). 結果
Tet10 の性能計測結果を図 2.13.に示す。CPU ベースの「富岳」では高いメモリバンド幅を活用した
BCRS 形式の行列ベクトル積が有効であり、「富岳」のアーキテクチャを考慮して構築していない一般的な BCRS 実装である BCRS-OpenMP でも演算ピークの 2.9%と FEM の行列ベクトル積カーネルとして 一 般 的 な 性 能 が で てお り 、 さ ら に 、 SIMD デ ー タ ア ク セ ス や プ リフ ェ ッ チ を 適 切 に 考 慮し た
BCRS(Fugaku)-OpenMP により演算ピーク性能比 6.5%とメモリバンド幅の 57%という高い性能が得られている。その一方で「富岳」においては「富岳」用にチューニングされた EBE を使っても BCRS より- 68 -

も遅く、メモリ使用量削減や行列構築コストの削減といった EBE 法の利点はあるものの、行列ベクトル
積単体における計算速度での恩恵はない結果となっている。その一方で、CPU であっても近年のコア性能の向上による相対的なメモリバンド幅比の低下が顕著な Xeon CPU では傾向が異なり、ABCI の CPUを使った場合、メモリバンド幅を活用できている BCRS-OpenMP, BCRS-OpenACC（それぞれ、メモリ
バ ン ド 幅 の 55.4%, 74.0% ） よ り も SIMD/ マ ル チ コ ア に 適 し た EBE (EBE(Fugaku)-OpenMP,
EBE(Fugaku)-OpenACC) を使った方が速いという結果になった。これは、EBE によるメモリ転送量を演算量に転嫁するアルゴリズムが有効に働いているためと考えられる。ただし、同じ EBE アルゴリズム
でも、高速な atomic add 性能を前提に作られた EBE(GPU)-OpenACC は BCRS よりも遅くなっており、SIMD/マルチコアを適切に活用しない方法を使うと同じ EBE アルゴリズムでも性能が得られない
ことがわかる。このように、CPU アーキテクチャにおいては従来からの BCRS-OpenMP による方法が有効であり、一部のアーキテクチャにおいては実装が少し複雑になるものの EBE により高速化される場合もあることがわかる。

図 2.13 Tet10 での測定結果
括弧内に FP32 FLOPS、メモリバンド幅のピーク比を示す。

- 69 -

図 2.14. Tet4 での測定結果
括弧内に FP32 FLOPS、メモリバンド幅のピーク比を示す。

GPU においてはメモリバンド幅・ピーク演算性能ともに CPU より大きく、適切な実装により
BCRS/EBE ともに大幅な高速化が期待される。なお、例えば「富岳」１ノードと A100 GPU １枚を比較した場合で、メモリバンド幅は 1.5 倍に対し演算ピークは 3.2 倍と、演算ピークの向上幅が大きく、特に EBE の有効性が高まると期待される。実際、ABCI 上で GPU を使った場合、BCRS-OpenACC は
CPU のみを使った場合と比較して、CPU と GPU のメモリバンド幅比の 1:30 に近い約 22 倍の高速化となっている(バンド幅性能から少し下がっている要因として疎行列の非零成分に対するアクセスが
coalesced アクセスではないことに起因する可能性がある)。
EBE(GPU)-OpenACC を使用することで、A100 GPU 上での BCRS-OpenACC のさらに 3.4 倍の高速化が実現しており、CPU 上で BCRS-OpenACC を用いた場合の 77 倍の高速化となる。このように、行列ベクトル積カーネルにおいて GPU を使うことによって、BCRS/EBE ともに高速化できることを示しており、とくに GPU アーキテクチャに適した EBE を使うことで大幅な高速化が期待できることがわかる。なお、今回の BCRS の実装が適切であることを確かめるため、cuSPARSE ライブラリの疎行列ベクトル積関数を使った BCRS-cuSPARSE (GPU 上で実行)と手製の BCRS-OpenACC (GPU 上で実行)の性能を比較したが、BCRS-cuSPARSE の性能の方が低く、手製の BCRS 実装で妥当な性能が得られていることを確認している（cuSPARSE 関数においては y ⇐ αAx + βy を計算しており、行列ベクトル積
以外の計算も行っていることに起因していると思われる。今回の計測ではα=1、 β= 0 とすることで行
列ベクトル積を計算している）。また、メモリ格納型の行列ベクトル積アルゴリズムにおいて本研究ではBCRS を用いたが、一般的に BCRS 以外の格納形式はメモリアクセス量を多少増やしても演算スケジュ
ーリングやメモリアクセスパターンの改善を目的として作られており、BCRS よりもメモリアクセス量は増加することとなる。今回の対象疎行列(Tet10)では、非零成分数は約 630 million/GPU であり、BCRS
形式で A を保持して行列ベクトル積を計算した場合のメモリアクセス量の下限値は 2.8 GB 程度である。
A100 GPU のメモリバンド幅は 1555 GB/s であることから、仮にこのメモリバンド幅を 100%使えたとしても、このメモリアクセスには最低でも 1.8 ms は必要である。つまり、メモリ格納型の行列ベクトル- 70 -

積アルゴリズムでは、1.8 ms 以上は必須となる。一方で、EBE 法を使った場合(EBE(GPU)-OpenACC
on GPU)の計算時間は 1.13 ms であり、疎行列をメモリに格納する場合、EBE より速くなることはない。
このように物理シミュレーションにおいて GPU、特に EBE を使用することで高速な解析が実現できることが期待される。
最後に、要素次数を 1 次に下げた Tet4 での各種カーネルの性能を計測した。図 2.14.に示すように、
このカーネルにおいても GPU は BCRS, EBE ともに CPU より高速化しており、Xeon CPU 比でそれぞれ 21 倍と 10 倍の高速化が得られている。Tet4 の場合では Tet10 の場合と異なり、「富岳」
・Xeon CPU・
A100 GPU のどのアーキテクチャにおいても、Tet10 での計測結果と比べて BCRS の優位性が増し、「富岳」では BCRS が EBE の約 3 倍、GPU においても BCRS が EBE の約 2 倍高速という結果となった。
これは、③-(d)-2). で想定していた通りの結果といえる。なお、ABCI (CPU)は「富岳」
・GPU と比べてハードウェアの演算ピークに対するメモリバンド幅比が相対的に低いために EBE と BCRS がほぼ同等の性能となっている。このように、要素タイプによる演算密度の変化により適切な行列ベクトル積アルゴ
リズムも変化することがわかるが、主に低次要素で BCRS、高次要素で EBE と使い分けることで幅広い問題で GPU を高速に活用できると期待される。
③-(d). 適用例
本研究で得られた知見を、液状化を考慮した地震応答解析に利用し、GPU を使用した解析の高速化を実現した。アプリケーションの開発にあたっては、シングルコードで開発を行うことが望ましい。それには大きく 2 つの理由がある。1 つ目は、アプリケーションの機能の改良である。複数アーキテクチャについて別々のコードがある場合、同じような開発を複数行う必要があり、さらにはバグの温床にもなりえる。2 つ目は、計算機のアーキテクチャが今後も変わり続けていくことが予想されることである。新しいアーキテクチャが出てくるたびに、新たにコード開発を行うのは現実的ではない。
本研究では、OpenACC を使った実装により、GPU によるシミュレーションの高速化を実現した。GPUを使って解析を高速化するには、GPU 用のコード開発が必要である。GPU コードの開発には、CUDAが良く使われてきたが、高開発コストであり、CPU 上で実行できなくなるというデメリットもある(デバッグなどは CPU で行う方が容易なため、CPU 上で実行できないと、コードの開発コストは格段に増加する)。OpenACC を使う場合、CPU 用のコードに指示文を挿入するだけで GPU 用のコードが開発できるため、CUDA を使う場合よりも開発コストが低い。さらに、コンパイラが対応していれば、コンパイ
ルオプションを変えるだけで GPU 用にもマルチコア CPU 用にもコンパイル可能である。そのため、GPU で高速化を実現しつつ、CPU でもスレッド並列を使用した計算が可能である。また、仮にコンパイ
ラが対応していなくても、シングルコア用のコードとして正常に動作する。OpenACC を使用すること
で、CPU/GPU の双方で動くパフォーマンス・ポータブルなシングルコードのアプリケーションを実装することができた。
我々は、液状化を考慮した地震応答解析の手法として、Intel Xeon Phi CPU 向け、GPU 向けの開発をおこなってきた。しかし、これらはそれぞれの対象アーキテクチャを主眼としていたため、別のアーキ
テクチャでは性能が出なかったり、動かなかったりする。本研究における開発により、これらを別々に管理する必要がなくなり、より効率的なアプリケーションの使用・改良が可能になると期待される。

- 71 -

③-(d)-1). 対象問題
対象問題は液状化を考慮した地震応答解析である。地盤の液状化とは、地震時に地盤がせん断剛性を失い、液体のように振舞う現象である。地盤の液状化は、建物の傾きの原因となる地盤の不同沈下や、マンホールなどの地下構造物の浮き上がり、噴砂などの被害を引き起こす。それらの被害を推定し、対策工の性能を評価するには、液状化を考慮した地震応答解析が必要であり、その高速化は非常に重要である。
地震応答解析では、式(1)を 104−5 回解くこととなる。液状化を考慮した解析においては、要素係数行列
Ae は以下のように計算される。

ただし、α, βは要素ごとに決まるスカラー定数である。また、M は要素質量行列、B は節点変位をひずみに変換するための行列であり、要素の形状によって決まる。行列 D は地盤の弾塑性行列であり、6×6の行列である。本研究では、液状化を考慮するため、Iai et al. (1992)の地盤構成則に従って、弾塑性行列を以下のように計算している。

ここで、Ku と Kf は地盤と間隙水の体積弾性係数、n は間隙率である。γi 、 wsi 、τi 、ni は 1 次元ばねのひずみ、重み係数、1 次元ばねの応力、3 次元ひずみを 1 次元ひずみに変換するためのベクトルである。Ns は 1 次元ばねの本数であり、本研究では 300 である。また、m = {1 1 1 0 0 0}T である。ここで、dτi/dγi は液状化の進行を考慮した経験式から計算される。
③-(d)-2). ソルバー式(1)のソルバーとして、本研究では、可変的前処理付き共役勾配法に基づく手法を用いた。可変的前処理付き共役勾配法は前処理付き共役勾配(PCG) 法の一種で、各反復で、「方程式 Az = r を粗く解く」
ことを前処理としている。前処理での Az = r の求解にはブロックヤコビ前処理を用いた共役勾配法を使用する。可変的前処理付き共役勾配法では、前処理が求解の計算コストの多くを占める。前処理において、マルチグリッド法を使用することで計算コストを削減している。求解対象となるモデルは、四面体 2次要素からなる地盤構造物モデルであるが、前処理では、Az = r を四面体 1 次要素のモデルにマッピン
グした Aczc = rc をまず解き、その解を四面体 2 次要素のモデルにマッピングしたものを初期解として、
Az = r を解く。低次モードを計算コストの低い 1 次要素モデル上で求めることで、計算コストの高い 2
次要素モデル上での反復回数を削減でき、計算コストが削減されている。前処理で、Aczc = rc を解くための反復を inner coarse loop、Az = r を解くためのループを inner fine loop と呼ぶ。また、これらの前処
理を各反復で行うループを outer loop と呼ぶ。前処理の計算コストをさらに削減するために、前処理での計算はすべて FP32 変数で行う(それ以外の計算は FP64 で行う)。近年の計算機では、FP64 変数に対する演算よりも FP32 変数に対する演算の方が高速であることが多いため、計算の高速化が期待できる(本研究で使用する A100 GPU では、peak FP32 FLOPS は peak FP64 FLOPS の 2 倍である)。また、- 72 -

FP32 を使う場合、メモリアクセス量や MPI 通信における通信量も半分となるため、計算時間の短縮につながると期待される。前処理は outer loop の反復回数を減らすためのものであり、最終的な解の精度には影響しないため、FP32 を前処理で使ったとしても、解の精度は変わらない。以上のアルゴリズムを
Algorithm 1 に示す。

③-(d)-3). 実装
行列ベクトル積に BCRS で疎行列を格納するアルゴリズムを採用し、CPU 上で解析を行うコード(CPU-BCRS) をベースラインとし、それを GPU ポーティングしたコード(GPU-BCRS)を開発した。
CPU-GPU 間のデータ転送を最小限にするため、基本的にすべての計算を GPU 上で行い、MPI 通信を行うデータと解析結果として出力するデータのみを GPU-CPU 間で転送する実装とした。さらに、四面体 2 次要素モデル上での行列ベクトル積では、EBE 法が GPU で効果的であることが③―(c)で示されたため、inner fine loop と outer loop における行列ベクトル積に EBE 法を用いたコード(GPU-EBE)を開発して、さらなる高速化を目指した。CPU-BCRS と GPU-BCRS はほとんど同じコードであり、GPU-EBE も行列ベクトル積以外は共通である。OpenACC を使用することで、CPU/GPU の双方で動くパフォーマンス・ポータブルなシングルコードのアプリケーションを開発することができた。

- 73 -

③-(d)-4). 液状化を考慮した地震応答解析における行列ベクトル積の計算
液状化を考慮した地震応答解析の場合、要素剛性行列 Ae の計算が通常の(液状化を考慮しない) 地震
応答解析とは若干異なる。そのため、ABCI の 1 計算ノード(A)上で追加の性能計測を行った。問題設定
は、③-(c)と同じである。図 2.15.に結果を示す。液状化を考慮した場合でも、EBE の方が BCRS で疎行列を格納する場合よりも 2.4 倍速いことが確認された。また、ピーク FP32 FLOPS の 17.6 %という高い性能を達成している。

図 2.15. 四面体二次要素における液状化解析（EBE-LIQ(GPU)-OpenACC）用のEBE 行列ベクトル積の性能計測結果
BCRS-OpenACC, EBE(GPU)-OpenACC は図 2.13.の再掲。カッコ内の数字はそれぞれ FP32 peak 及びメモリバンド幅ピーク比を示す。

③-(d)-5). 構成則の計算の性能計測
液状化を考慮した地震応答解析では、構成則の計算は解析全体の 20％程度の計算時間を占める。構成則の計算では、各要素の弾塑性行列の計算をしているが、その際のメモリアクセスパターンが 2 通り考
えられる(図 2.16.)。

図 2.16. 弾塑性行列の計算におけるメモリアクセスパターン
配列は column major order (Fortran style)で記述。

- 74 -GPU と CPU でそれぞれのパターンで計測を行った。GPU ではパターン 1 がパターン 2 の 3.9 倍速
かった。これは、パターン 2 では parameters へのアクセスがストライドアクセスとなり、DRAM アクセスが増加したためと考えられる。GPU では、パターン 1 はピーク FP64 FLOPS の 37.3 %という高い性能を達成した。一方で、CPU ではパターン 2 がパターン 1 の 1.5 倍速かった。CPU と GPU で性能の高いメモリアクセスパターンが異なることが分かった。一方で、なるべくシングルコードでの開発を行う
ことが望まれる。そこで、2 次元配列を 1 次元化し、インデックスの計算部分においてのみ#ifdef を使うことで、CPU, GPU に適したメモリアクセスパターンとなるようにした（図 2.18.）。CPU-BCRS と GPUBCRS で異なるのはここのみで、それ以外は全く同じコードである。

図 2.17. 構成則の計算の性能計測

図 2.18. 弾塑性行列の計算におけるメモリアクセスパターンの切り替え
配列は column major order (Fortran style)で記述。

③-(d)-6). 解析全体の性能計測
性能計測用の地盤モデルを使用して、解析全体の性能計測を行った。モデル自由度は 58,574,715、要
素数は 14,434,560 である。ABCI 2 計算ノード(A) 上で 16 MPI プロセスで 100 時間ステップの液状化- 75 -

解析を行った。計測したケースは、③―(d)―3）で説明した CPU-BCRS, GPU-BCRS, GPU-EBE (いずれも MPI と OpenACC による並列化) に加え、スーパーコンピュータ Oakforest-PACS(OFP)用に MPIと OpenMP で並列化したコード(CPU-OFP) の 4 ケースである。CPU-BCRS, CPU-OFP は CPU のみを使用して実行し、GPU-BCRS, GPU-EBE は GPU を使用して実行した。
図 2.19.に性能計測結果を示す。まず、ベースラインである CPU-BCRS は CPU-OFP の 1.16 倍の計算時間だった。CPU-OFP は OFP (Intel Xeon Phi CPU を搭載)用に開発されたものであり、ABCI (Intel
Xeon CPU を搭載)用のチューニングにより、高速化できる可能性はある。しかし、特別なチューニング
を行っていない OpenACC コードである CPU-BCRS が既往の CPU 用コードとおおむね同程度の性能を出せることが示されており、マルチコア CPU 用の実装における OpenACC の有効性が示されたといえる。GPU-BCRS は CPU-BCRS と比較して 24 倍、CPU-OFP と比較しても 21 倍の高速化を実現した。
GPU-EBE は GPU-BCRS と比較して 1.7 倍速く、CPU-BCRS と比較すると 42 倍の高速化である。GPU
の特徴を考慮して、演算量が多いがメモリアクセス量は少ないアルゴリズムを使うことで GPU の性能をより有効に活用し、解析を高速化することができたといえる。

図 2.19. 液状化を考慮した地震応答解析の計算時間の比較

これらの結果は、OpenACC を使った GPU ポーティングにより、飛躍的な高速化ができたことを示している。また、GPU-BCRS のコードは、既往の CPU 手法と同程度の性能が出る CPU-BCRS と同じコードであり、OpenACC を使うことで、パフォーマンス・ポータブルなシングルコードでのアプリケーシ
ョンの開発が実現した。
③-(d)-7). 実地形を模した地盤モデルについての液状化解析
最後に、開発コード(GPU-EBE) を使用して、河口付近の地形を模した地盤モデルについての液状化解析を行った。図 2.20.にモデルの概観を示す。モデル自由度は 89,146,716, 要素数は 21,616,202 である。
入力波としては、平成 7 年の兵庫県南部地震における観測波を用いた。ABCI の 6 計算ノードを使用して、時間ステップ幅 dt=0.001s で 30,000 時間ステップの計算を行った。図 2.21. のような変位応答が得られた。

- 76 -

図 2.20. 実地形を模した地盤モデル

図 2.21. 解析結果（図 2.20.の点 A における変位応答）

計算時間は 1 時間 58 分であった。この問題を CPU で計算した場合、3.4 日かかることになる(CPU に
おいては GPU を用いた場合の 42 倍の時間がかかると仮定)。GPU を使うことで、格段に速く解析ができるようになった。地盤や入力地震動の不確実性を考慮するためには多数回の解析が必要である。1 回の解析に何日もかかってしまう場合、そのような多数回解析は困難であったが、GPU による高速化により多数回の解析が容易にできるようになると期待される。
③-(e). まとめ
本研究では「富岳」で開発した高効率行列ベクトル積手法の汎用化のための開発を実施した。ここでは有限要素法における疎行列ベクトル積を「富岳」だけでなく GPU 計算機においても高速実行できるよう、GPU に適した疎行列ベクトル積アルゴリズム・実装を開発した。1 ABCI 計算ノード上で性能計測を行ったところ、BCRS 形式で行列を保持し計算に用いる一般的な疎行列ベクトル積アルゴリズムの場合でGPU 版は CPU 版の 22 倍の高速化を達成した。また、GPU の高い演算性能をより有効活用できる EBE
法(演算量は多いがメモリアクセス量が少ない)を用いることによって、CPU の 77 倍高速化することが分かった。
この知見をもとに GPU を使って液状化解析を高速化した。OpenACC を用いて実装を行い、GPU を使うことで 42 倍高速化した。OpenACC を使うことで GPU による解析の高速化が実現可能であり、かつ、CPU と GPU に対してパフォーマンス・ポータブルなコード開発が可能なことが示された。

- 77 -

④ 地震観測データをシミュレーションにて有効活用するための高速データ解析手法の開発
本プロジェクトでは、与えられた地盤モデル等の地下構造に対して、地震動や地殻変動の大規模高速計
算を行うコード開発を主に行っているが、開発されたコードが実際の問題で有効に活用されるためには、
地盤モデルがより現実的であることが必要となる。地盤モデルの改良には高コストなデータ前処理が必要となるが、その計算は NVIDIA GPU の行列積演算加速ユニットである Tensor Core を用いることで高効率で実行できると見込まれる。データ前処理部を GPU で高速計算できるようになればその分の「富岳」計算資源を地震動や地殻変動の大規模高速計算に充てることができるようになるため、本研究においてはこのデータ処理の NVIDIA GPU における高速化を検討した。近年では多数点での観測により従来よりも格段に大容量の地震観測データが入手可能となっており、これを数値シミュレーション・データ分析することで有用な情報が引き出せると期待されている。例えば地震波干渉法では、異なる 2 点で観測された応答波形の相互相関関数を求め、長期間に渡りスタッキングすることで、観測点間の媒質情報(地盤情報)を抽出し地盤モデルの改良に用いることができるが、観測点が多数・波形が長期間にわたる場合は計算量が膨大になりデータ分析のボトルネックとなっていた。なお、相互相関関数は時系列データを扱う
多くの分野で登場する計算であり、近年多数のデータが蓄積されているという背景から、この計算の高速化は他の分野においても有用と期待される。相互相関関数計算は行列-行列積演算の形となるため、近年の NVIDIA GPU の行列積演算加速ユニットである Tensor Core を用いることで大幅な高速化が見込まれる。そこで本研究では Ampere アーキテクチャから使用可能となった TensorFloat-32 と呼ばれる新たな precision data type を使用した Tensor Core 演算を用い、相互相関関数計算と Tensor Core の両特性を考慮した高速な相互相関関数の計算手法を開発した(図 2.22.)。ここではまず Ampere アーキテクチャ
における TensorFloat-32 を用いた Tensor Core 演算を対象とし、shared メモリと fragment のマッピングを調べた。この調査結果をもとに PTX インラインアセンブリを用いることで通常の API を介さない
低レベル記述の実装を実現し、メモリ帯域に律速されにくい手法の計算手法を開発した。地震波干渉法における相互相関関数計算を対象とした性能計測において、本手法は A100 GPU 1 枚あたり 53.56 TFLOPSという非常に高い実行性能を得た(表 2.8.参照)。これは一般的に GPU 上で広く利用されている線形代数
演算ライブラリである cuBLAS を用いて TensorFloat-32 (Tensor Core)で計算した場合と比較し 5.97 倍
高速であった。64-bit 浮動小数点演算にて計算した結果と比べ、提案手法での計算結果の精度は許容範囲に十分収まっており、TensorFloat-32 を用いた Tensor Core 演算の干渉法への適用可能性が示された。
本研究での高速化検討により、地震観測データをはじめとした多数のデータをより有効に活用できるよ
うになると期待される。

- 78 -

図 2.22．Tensor Core を用いた TensorFloat-32 の相互相関関数計算の概要
ここでは、1 スレッドブロックの計算の概略を示している(L は相互相関関数の総時間シフト量を示す)。

表 2.8．各カーネルの計算性能と精度 (V100 16GB SXM 1 枚、A100 40GB SXM 1 枚、及び、Xeon
Platinum 8360Y (36 cores)を利用した結果)
Yamaguchi (2019)は前世代の GPU である V100 GPU 向けに最適化された相互相関関数計算手法の結果。

（３）開発した計算アプリケーション群を実務で活用するための環境整備（１）において整備している国の被害想定のためのハザード評価に使われるものと同等の長周期地震
動計算アプリケーション群の実務への活用を進めるため、連携機関の各企業が独自の問題設定で実施するためのサポート体制を、海洋研究開発機構が中心となって整備している。対象としているのは計算アプ
リケーション群の基幹数値解析コンポーネントである E-wave FEM の計算環境である。課題実施期間の三年間で、E-wave FEM の機能・ユーザビリティの向上と連携企業による活用のサポートを重点的に進めてきた。本稿①にて、本年度におけるこれらの取り組みについて述べる。また連携企業は、整備した Ewave FEM 計算環境を用い、企業でニーズのある地震動計算への適用を進めた。その中で本年度に得られた成果について、②で述べる。

- 79 -

① E-wave FEM の機能・ユーザビリティの向上と連携企業による活用のサポート①-(a) 実施した計算環境整備(i) 計算機システムの運用停止・開始への対応


概要：E-wave FEM の計算環境は、大規模共有メモリ計算機からなる計算環境（以下、LM システ
ムと記載 ）を計算モデル（メッシュ）生成・並列計算用領域分割用に、高並列スーパーコンピュータからなる計算環境（以下、SC システムと記載）を高並列地震動計算用に、連接して利用する。LM・SC システムともに、複数の計算機システムをサポートしており、計算機の運用停止・開始に合わせたアップデートを実施した。



従来対応していた計算機のサポート停止：令和４年度は、「富岳」と共に SC システムとして利用可
能としていた Oakforest-PACS (東京大学情報基盤センターと筑波大学計算科学研究センターが共同運営する、最先端共同 HPC 基盤施設 JCAHPC：Joint Center for Advanced High Performance
Computing の共同利用スーパーコンピュータ)が令和4年3月に運用停止となった。また、LM シス
テムを構成する大規模共有メモリ計算機 Jasper を海洋研究開発機構から理化学研究所へと移管した。それに対応して、Oakforest-PACS と Jasper に対するサポートを停止した。



計算機の新規導入：これらの計算機システムの運用停止に伴い、令和4年度は以下の計算機に対するサポートを新たに開始した。LM システムとして、第2階層スーパーコンピュータである地球シミ
ュレータ第4世代 (ES4:海洋研究開発機構: https://www.jamstec.go.jp/es/jp/system)を使用可能とした。ES4は VE（Vector Engine）ノード、CPU ノード、GPU ノードにより構成されるヘテロジニ
アスアーキテクチャであるが、このうち大規模メモリを有する GPU ノードを使用可能とした。また、海洋研究開発機構にて運用開始した大規模共有メモリ型サーバ（virgo）にも対応可能とした。
SC システムとしては、「富岳」に加えて、ES4の VE ノードおよび CPU ノードも使用可能とした。



現在使用可能な計算機システムのまとめ：E-wave FEM のプリプロセス（LM システム利用）において ES4（GPU ノード）
、virgo が、FEM 本解析とポストプロセス（SC システム利用）において「富岳」および ES4（VE および CPU ノード）が利用可能となった。なお、別々の計算機を利用する際、使用感が統一的となるように配慮している。

(ii) 計算環境に対する機能の付加、バグ修正


概要：令和4年度は、実地盤構造に即した現実的な設定における地震動計算に対するユーザビリテ
ィを高めるための機能を追加した。



実地盤データの読み込みに対応した計算環境：E-wave FEM による地震動計算では、公開されている実地盤構造データを読み込んで計算モデルを作成することが多い。日本国内においては、2種類のデータセット、すなわち全国1次地下構造モデル（暫定版：
https://www.jishin.go.jp/evaluation/seismic_hazard_map/lpshm/12_choshuki_dat)、および JSHIS(地震ハザードステーション：防災科学技術研究所：https://www.jshis.bosai.go.jp/dstrctmap-2020)が広く使われる。これらのデータを計算に用いる際は、地理座標系で定義される地盤層境界の形状データ座標値を、計算モデルの定義されるデカルト座標系におけ- 80 -

る座標値に変換したうえで、必要な領域を切り出す必要がある。これらの作業を自動で行う計算環
境を実装した。


E-wave FEM 複数回計算に対応した計算環境：E-wave FEM の用途の一つとして、同じ地盤デー
タに対して複数の震源パラメータを入力した多数回計算があることが分かった。そこで、複数組の震源パラメータ・観測点情報に対する地震動をまとめて設定・計算可能な計算環境を構築した。複数通りの計算に対して各々にサブディレクトリを割り当て、複数回計算選択のパラメータを設定することで、単独計算の場合と同じ計算ステップ、同じインターフェースで、複数回計算を実現できる機能を導入した。



ユーザビリティのさらなる改善：多様な計算機システムを利用対象としたことで、ユーザに求められる並列計算のパラメータ設定は、より複雑になっている。そこで、適切なパラメータ設定（FEM
本解析における MPI プロセス数）をログに表示し、計算以前にユーザに提案する機能を追加した。具体的には、SC システムにおける FEM 本解析（MPI 計算）では、適切な MPI プロセス数に分割されたメッシュモデルを利用する。LM システムで実施されるメッシュモデル分割の前に、メッシュのデータサイズに応じて適切な MPI プロセス数（ノード数）をログに記載しユーザに提案する機能を付加した。また、計算規模拡大のため出力される中間ファイル（特に、LM システムで生成される分割前のメッシュモデル）のサイズも大型化している。SC システムではこの分割前のメッシュモデルも必要となるが、その大型化のため LM システムから SC システムへの転送が不可となる、或いは時間がかかることが判明した。そこで、両システム間でファイル転送するデータ量を削減するために、元々SC システム上で分割前メッシュモデルを必要としていた作業を LM シス
テムでも実行できるような拡張機能を追加した。



細かなバグ修正に伴う8回のアップデート：これらの様々な機能追加だけでなく、細かいバグ修正を含めて令和4年度は、8回のアップデートを行なった。また、2回の主要アップデートの段階でユーザへの連絡を行なった。

(iii) E-wave FEM 計算モデル断面データの抽出・可視化のための計算環境の構築


概要：連携企業からの参加者は、既に他の地震動計算手法の使用経験を有することがほとんどである。そのため、E-wave FEM の企業での利用においては、まずそれらの既存のツールとの結果比較が行われれることが多い。その一環として、両ツールで使用する計算モデルを任意の二次元断面において比較するニーズがあることが判明した。そこで、予め定めた E-wave FEM の計算領域内において、計算に導入する元となる実地盤データ、実地盤データから計算設定に応じて切り出した、
地盤入力データ、地盤入力データに対し E-wave FEM 計算用メッシュを適用したメッシュモデ
ル、の3種類のデータおよびモデル（まとめて以下、3種データ）を抽出、可視化できる機能を開発し、実装した。本可視化システムを用いて作成した関東地方のP波速度構造断面図の例が、本報告書（１）②-(c)にて掲載する図1.2.28、1.2.29に示されている。



可視化断面と可視化要素の選択：本可視化システムは、断面に対する視線方向ベクトルを設定可能とすることで、任意の視線方向から任意の傾斜角の断面に対して3種データを抽出、可視化できる。デフォルト設定では、鉛直断面を測線と直行する水平方向から可視化する設定とした。また、
断面の縦横比の変更、断面の一部切り出し、拡大縮小表示を可能とした。更に、3種データから少なくとも2つを選択して、重畳表示あるいは単独表示を可能とした。重畳表示では、重ね合わせる- 81 -

順序も交換可能とした。また、可視化する2点間の組み合わせ数は、入力パラメータとして設定可能とした。表示色は、層内の物性、層境界線、メッシュに対して自由に設定可能であり、単独色（無色も含む）、カラーマップを使った連続色のいずれも設定できる。可視化要素として、予め各層に定めた物性値（p 波速度、s 波速度、密度、Q 値）から選択できる。


複数の可視化手順のサポート：本可視化システムは、LM システムにおいて実施する E-wave FEM
計算のプリプロセス（実地盤データの取り込み、実地盤データから地盤層データの抽出、地盤層データからメッシュモデルの作成終了後）の中で実施できる。同時に、全計算終了後の過去の保存データからの可視化も可能とした。その際は、補助スクリプトによりファイルリンクを自動化し、前記プリプロセス内で実施した場合と使用感を統一し、同じステップ、同じインターフェースで操作可能である。

①-(b) 連携企業のサポート
令和3年度までには、上記のような計算環境の整備・アップデートに加え、全6回のオンライン形式のチュートリアル、詳細なマニュアルの整備・改定等を実施してきた。利用が進んできた令和3・4年度には、海洋研究開発機構と連携企業ユーザとの間での研究の進捗に応じた定期的な個別打合せにも力を入れてきた。特に、トラブル発生時には海洋研究開発機構側でファイル設定を精査し、回避策の提案を行なった。また、計算システムの利用想定を超える用途に対しては、実現のための方策を提案した。
この様な個別打合せの中で明らかになった問題への回避策、解決策の提案は、計算システムの改良へとフィードバックされた。例えば、プリプロセスでの実地盤データの読み込み機能、および地震動計算における複数回計算機能の追加は、連携企業ユーザ側での試用（より現実的な実地盤モデルの導入、震源過程の不確定性を反映した多数回計算設定等）を通じた要望に端を発している。また、プリプロセスでの大規模データ転送機能および計算資源提案機能は、連携企業ユーザからの問い合わせに端を発している。計算モデルの断面可視化については、連携企業ユーザが比較対象となる別の数値計算手法に対して個別に実施していることもあり、本計算環境内での実現を強く要望されている。今回追加したこれらの機能は、連携企業ユーザから期待される機能であり、更に連携企業ユーザからのフィードバックを受けて、利便性を高めつつ機能強化を続ける予定である。
② 企業においてニーズのある地震動計算への適用：多様な震源シナリオを考慮した長周期地震動評価と建築物応答解析－日本近海のプレート境界地震を対象とした検討－
チュートリアル・本課題でのサポートに基づき、連携企業メンバーによって、企業でニーズのある地震
動計算への E-wave FEM の活用が進められている。
その中で本年度得られた成果の一つについて述べる。
地震災害に対してレジリエントな社会の実現に向けた対策の 1 つとして、被害拡大を防ぎ早期の復旧復興に寄与する耐震性能の高い建築物を設計することが挙げられる。現在、高い耐震性能が要求される超高層建築物では、建設地点の地震環境を考慮した設計用地震動、いわゆるサイト波を入力とする時刻歴応
答解析が行われている。これらの建築物は 1 次固有周期が長いため、サイト波の長周期成分（長周期地震動）の影響が大きく、適切な地震動評価が求められている。
文部科学省地震調査研究推進本部（以下、地震本部）では、活断層や海溝型地震を対象に、地震規模や- 82 -

発生確率を予測・公表する長期評価を行ってきた。従来までは、各震源域で同程度の規模の地震が繰り返し発生する固有地震の考え方に基づいていた。しかし、平成 23 年 3 月 11 日東北地方太平洋沖地震のような超巨大地震を想定できなかった背景を受け、これまでの長期評価の方法を見直し、最大クラスも含めた地震の多様性や不確実性を取り入れた評価が重要視され、震源特性の不確実性を考慮した多数の破壊シナリオの地震動評価が行われるようになった[1]。
固有地震の考え方に基づき震源特性や破壊シナリオを一意に決めた評価に対し、震源特性の不確実性を考慮した地震動評価は、その選択により予測結果に大きなばらつきが生じる。例えば、五十嵐ら[2]は震源特性の不確実性を考慮するモンテカルロシミュレーションにより、多数の地震動を作成し、ばらつきの大きさと建物に及ぼす影響を評価している。また、ばらつきをもつ多数の地震動をサイト波へ適用した例として、大震研委員会[3]では、複数の破壊シナリオを考慮した多数の地震動評価結果から、建築物の要求耐震性能に応じてサイト波に対して考慮するばらつきの大きさを選択する方法を提案している。
平成 15 年十勝沖地震では、震源から遠く離れた苫小牧において長周期地震動の被害が顕在化した。これは勇払平野の厚い堆積層の影響であり、長周期地震動が震源から被害箇所に至る 3 次元地盤の影響を強く受けることが改めて認識された。3 次元地盤の不整形性の影響を考慮して長周期地震動を適切に評価するには、3 次元有限差分法（3D-FDM）や 3 次元有限要素法（3D-FEM）のような 3 次元数値解析手法が必要になる。計算コストの観点から、これらの手法では不確実性を考慮した多数ケースの検討を現実的な計算時間で評価することが難しかった。しかし近年の計算機性能の向上により、3 次元数値解析による多数ケースの長周期地震動の検討例が出てきている。例えば、地震本部の長周期地震動予測地図 2016 年
試作版[4]では、相模トラフの地震を対象に多数のシナリオを想定した評価と統計値の算出が行われている。また本課題でも[5]も、熊本地震前震を対象に震源特性の不確実性を考慮した長周期地震動を評価し、地震動のばらつきを算出している。今後は長周期地震動についても、震源特性や破壊シナリオの不確実性を考慮した多数ケースの地震動とばらつきを評価し、大震研委員会[3]のように、建築物に必要な耐震性能に応じて考慮するばらつきの大きさを選択し、サイト波として用いることが望ましい。そのためには、
震源特性の不確実性が地震動や建物応答に与える影響、その統計値を把握しておくことが重要である。そこで本検討では、日本海溝、南海トラフの地震について、複数の震源特性、破壊シナリオを考慮した地震動評価を「富岳」を用いて行い、震源特性が地震動のばらつきに与える影響を分析した。
本検討では、日本海溝、および南海トラフの地震を対象とした。なお日本海溝は、岩手県沖～青森県沖と、茨城県沖～岩手県沖（3.11 タイプ）の 2 震源域を対象とした。なお、地震動評価を昨年度実施した
千島海溝[6]に関しても参考に記載する。
地震動評価には、地震動計算用大規模 3D-FEM 計算プログラム「E-wave FEM」[7][8][9]を用いた。
地盤モデルは深部地盤に J-SHIS-V3.2[10]を用い、地震本部の全国 1 次地下構造モデル[11]のプレート構
造を下部に接続した。地盤モデルの範囲を図 3.2.1 に、物性を表 3.2.1 に示す。対象周期 3.3 秒以上、1 波
長 5 要素以上となるよう要素分割を行った。節点数、要素数および 1 ケースあたりの「富岳」ノード時間を表 3.2.2 に示す。なお、要素分割には ES4 を用いた。
震源位置は内閣府[12]や、川辺(平成 25 年) [13]を参考に決定した。震源特性は、震源域毎にレシピ[14]のスケーリング則に従うカスケードモデルとして設定した。また、表 3.2.3 中に示す震源特性に不確実性を考慮するため複数のケースを設定し、震源モデルを作成した。震源特性を表 3.2.4～3.2.7 に、震源モ- 83 -

デルを図 3.2.2～3.2.5 に示す。なお、走向角、傾斜角は太平洋プレート上面に沿うように設定し、すべり方向は太平洋プレートの進行方向とした。すべり速度時間関数は山本・吉村[15]と同様に図 3.2.6 に示す
2 つの三角形で表現した。
地震動評価には「富岳」の 1,024 ノード（49,152 並列）を用いた。また本検討では、地震動強さを擬
似速度応答スペクトル（pSv，h=5％）で評価することとする。
日本海溝（青森県沖～岩手県沖）の地震の全ケース平均 pSv を図 3.2.7 に、pSv のばらつき（自然対数標準偏差）を図 3.2.8 に示す。日本海溝（岩手県沖～茨城県沖）の地震の全ケース平均 pSv を図 3.2.9に、pSv のばらつきを図 3.2.10 に示す。南海トラフの地震の全ケース平均 pSv を図 3.2.11 に、pSv のばらつきを図 3.2.12 に示す。各震源域とも勇払平野、十勝平野、関東平野、濃尾平野、大阪平野など、堆積層が厚い地域で pSv が大きい。また、対象周期が長くなるほど、pSv が大きくなる傾向がみられる。
一方で、ばらつきの大きさは pSv の大小と必ずしも対応していないため、震源特性の違いが影響していると考えられる。
日本海溝（岩手県沖～茨城県沖）について、東北地方太平洋沖地震の川辺ほか（平成 25 年）[13]の震源モデルによる長周期地震動評価結果を、図 3.2.13 に示す。また、東北地方太平洋沖地震の K-NET、KiK-net の観測記録を図中に〇で示す。周期 5～7 秒に関しては、解析結果と観測記録はおおむね対応している。一方で、周期 3 秒では解析結果が観測記録より小さく、各点とも過小評価となり、周期 10 秒付
近では、特に関東平野の観測点で解析結果が観測記録より大きく、過大評価となっている。
周期 3 秒における過少評価は、すべり速度時間関数などの震源特性の設定や、断層破壊のすべり不均質を導入していない点等が影響していると考えられる。また周期 10 秒における過大評価は、アスペリテ
ィモデルを仮定した震源設定で評価できる周期として上限に近く、大すべり域を含む震源モデルによる地震動評価が必要となる可能性がある。本報では、破壊開始点、アスペリティサイズ、アスペリティ位置の 3 パラメータに不確実性を付与しているが、他のパラメータ設定に関する不確実性に関する検討は、
今後の課題である。

- 84 -

図 3.2.1 モデル化の範囲- 85 -

表 3.2.1 地盤物性
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

Vp
1600
1600
1700
1800
1800
2000
2000
2100
2100
2200
2300
2400
2400
2500
2500
2600
2700
3000
3200
3400
3500
3600
3700
3800
4000
4000
5000
4600
5000
5500
5800
6400
7500
5000
6800
8000
5400
6500
8100

Vs
350
400
450
500
550
600
650
700
750
800
850
900
950
1000
1100
1200
1300
1400
1500
1600
1700
1800
1900
2000
2100
2100
2700
2900
2700
3100
3400
3800
4500
2900
4000
4700
2800
3500
4600

ρ
1850
1850
1900
1900
1900
1900
1950
2000
2000
2000
2050
2050
2100
2100
2150
2150
2200
2250
2250
2300
2300
2350
2350
2400
2400
2400
2500
2550
2500
2600
2700
2800
3200
2400
2900
3200
2600
2800
3400

Q値
60
60
60
60
60
100
100
100
100
100
100
100
100
150
150
150
150
150
150
150
150
150
150
200
200
200
200
200
200
300
400
400
500
200
300
500
200
300
500

備考

J-SHIS 深部地盤モデル

全国 1 次地下構造モデル

上部地殻第 2 層
下部地殻
マントル
海洋性地殻第 2 層（フィリピン海プレート）
海洋性地殻第 3 層（フィリピン海プレート）
海洋性マントル（フィリピン海プレート）
海洋性地殻第 2 層（太平洋プレート）
海洋性地殻第 3 層（太平洋プレート）
海洋性マントル（太平洋プレート）

表 3.2.2 節点数・要素数および 1 ケースあたりの「富岳」ノード時間
モデルサイズ

節点数

要素数

1 ケースあたりの
ノード時間

日本海溝
（岩手県沖～青森県沖）630.4km×480.0km×70.0km
最小メッシュサイズ 200m

22.3 億16.5 億約 2000NH

日本海溝（3.11 タイプ）

700.8km×460.8km×70.0km
最小メッシュサイズ 200m
1049.6km×480.0km×60.0km最小メッシュサイズ 200m

23.2 億17.2 億

約 2100NH

34.6 億25.7 億

約 4400NH

南海トラフ- 86 -

表 3.2.3
破壊開始点
日本海溝（岩手県沖～青森県沖）
日本海溝（3.11 タイプ）

5 か所

南海トラフ

5 か所

震源域名

岩手県沖
～
青森県沖

1
2
3
背景領域
合計

面積
[km2]
2775
3318
1706
62425
70225

震源域名

岩手県沖
～
宮城県沖

福島県沖
～
茨城県沖

1
2
3
背景領域
合計
1
2
背景領域
合計

面積
[km2]
2425
1481
388
52900
57194
706
850
32131
33688

震源域名

日向灘

南海 1

南海 2

東南海 1

東南海 2

東海

1
2
背景領域
合計
1
2
背景領域
合計
1
2
背景領域
合計
1
2
背景領域
合計
1
2
背景領域
合計
1
2
背景領域
合計

面積
[km2]
913
1044
20419
22375
1469
1469
38894
41831
1594
825
25444
27863
825
1175
19944
21944
1044
688
17988
19719
481
456
7056
7994

20 ケース24 ケース30 ケース
5 ケース
28 ケース

Ar

面積

2.2
2.4
1.6
3.8

[km2]
5494
5856
3338
55538
70225

アスペリティが大きいモデル
平均
ライズ
地震
モーメント すべり量
タイム
[Nm]
[m]
[s]
3.81E+21
16.9
13.55
4.19E+21
17.5
13.99
1.80E+21
13.2
10.56
1.36E+22
6.0
24.96
2.34E+22
8.1

Ar
2.8
2.7
2.4
3.8

震源特性（日本海溝・3.11 タイプ）

アスペリティが小さいモデル
地震
平均
ライズ
モーメント すべり量
タイム
[Nm]
[m]
[s]
1.68E+21
16.9
9.00
8.01E+20
13.2
7.03
1.07E+20
6.7
3.60
1.46E+22
6.7
34.84
1.72E+22
7.3
3.10E+20
10.7
4.86
4.09E+20
11.7
5.33
7.06E+21
5.4
38.48
7.78E+21
5.6

表 3.2.6
アスペリティ
番号

浅い・深い
2 ケース

ケース数

震源特性（日本海溝・岩手県沖～青森県沖）

アスペリティが小さいモデル
平均
ライズ
地震
モーメント すべり量
タイム
[Nm]
[m]
[s]
1.87E+21
16.4
9.63
2.44E+21
17.9
10.53
8.99E+20
12.9
7.55
1.82E+22
7.1
24.96
2.34E+228.1

表 3.2.5
アスペリティ
番号

大きい・小さい
2 ケース

7 か所

表 3.2.4
アスペリティ
番号

アスペリティサイズ
アスペリティ位置
大きい・小さい
浅い・深い
2 ケース
2 ケース
大きい・小さい
浅い・深い
2 ケース
2 ケース
大きい・小さい
浅い・中位・深い
2 ケース
3 ケース
内閣府被害想定の配置 1 ケース

6 か所

千島海溝（昨年度検討）

考慮する震源特性の不確実性

Ar
2.2
2.4
1.6
3.8
2.2
2.4
3.8

面積
[km2]
4831
2906
763
48694
57194
1406
1625
30656
33688

アスペリティが大きいモデル
地震
平均
ライズ
モーメント すべり量
タイム
[Nm]
[m]
[s]
3.35E+21
16.9
12.70
1.56E+21
13.1
9.85
2.10E+20
6.7
5.05
1.21E+22
6.1
34.84
1.72E+22
7.3
6.25E+20
10.8
6.85
7.76E+20
11.6
7.37
6.38E+21
5.1
38.48
7.78E+21
5.6

Ar
2.8
2.7
2.4
3.8
2.8
2.7
3.8

震源特性（南海トラフ）

アスペリティが小さいモデル
地震
平均
ライズ
モーメント すべり量
タイム
[Nm]
[m]
[s]
3.31E+20
8.6
6.04
4.05E+20
9.2
6.46
3.48E+21
4.1
31.96
4.21E+21
4.5
7.56E+20
12.3
7.67
7.56E+20
12.3
7.67
9.26E+21
5.7
41.83
1.08E+22
6.1
7.41E+20
11.1
7.98
2.76E+20
8.0
5.74
4.84E+21
4.5
37.15
5.86E+21
5.0
2.76E+20
8.0
5.74
4.70E+20
9.5
6.86
3.35E+21
4.0
31.35
4.09E+21
4.4
3.99E+20
9.1
6.46
2.13E+20
7.4
5.24
2.87E+21
3.8
35.85
3.49E+21
4.2
1.10E+20
5.4
4.39
1.01E+20
5.3
4.27
6.89E+20
2.3
19.98
9.00E+20
2.7- 87 -

Ar
1.2
1.3
4.6
1.6
1.6
5.0
1.8
1.1
4.9
1.1
1.5
4.6
1.3
1.0
4.9
0.9
0.9
3.2

面積
[km2]
1788
2056
18531
22375
2881
2869
36081
41831
3156
1719
22988
27863
1625
2319
18000
21944
2019
1338
16363
19719
931
894
6169
7994

アスペリティが大きいモデル
地震
平均
ライズ
モーメント すべり量
タイム
[Nm]
[m]
[s]
6.48E+20
8.6
8.46
8.00E+20
9.3
9.07
2.77E+21
3.6
31.96
4.21E+21
4.5
1.49E+21
12.3
10.74
1.48E+21
12.2
10.71
7.81E+21
5.2
41.83
1.08E+22
6.1
1.46E+21
11.0
11.24
5.87E+20
8.1
8.29
3.81E+21
3.9
37.15
5.86E+21
5.0
5.44E+20
8.0
8.06
9.27E+20
9.5
9.63
2.62E+21
3.5
31.35
4.09E+21
4.4
7.71E+20
9.1
8.99
4.16E+20
7.4
7.31
2.30E+21
3.3
35.85
3.49E+21
4.2
2.12E+20
5.4
6.10
1.99E+20
5.3
5.98
4.89E+20
1.9
19.98
9.00E+20
2.7

Ar
1.9
2.1
4.6
2.5
2.4
5.0
2.5
1.9
4.9
1.8
2.2
4.6
2.1
1.5
4.9
1.2
1.2
3.2

表 3.2.7
震源域名

色丹島沖
及び
択捉島沖

根室沖

十勝沖

アスペリティ
番号
1
2
3
4
背景領域
合計
1
2
背景領域
合計
1
2
3
背景領域
合計

面積
[km2]
2413
2050
763
775
67819
73819
2131
775
20413
23319
700
1919
781
32969
36369

震源特性（千島海溝[6]・昨年度実施）

アスペリティが小さいモデル
平均
ライズ
地震
モーメント すべり量
タイム
[Nm]
[m]
[s]
1.92E+21
19.4
8.98
1.50E+21
17.8
8.27
3.40E+20
10.9
5.05
3.49E+20
11.0
5.09
2.11E+22
7.6
26.24
2.52E+22
8.3
9.16E+20
10.5
8.44
2.01E+20
6.3
5.09
3.37E+21
4.0
31.57
4.48E+21
4.7
2.43E+20
8.5
4.84
1.10E+21
14.0
8.01
2.87E+20
8.9
5.11
7.10E+21
5.3
37.87
8.73E+21
5.9- 88 -

Ar
2.1
1.9
1.0
1.0
4.0
1.9
1.0
4.6
1.0
1.8
1.0
5.0

面積
[km2]
4969
4188
1569
1575
61519
73819
4369
1569
17381
23319
1469
3975
1613
29313
36369

アスペリティが大きいモデル
平均
ライズ
地震
モーメント すべり量
タイム
[Nm]
[m]
[s]
3.95E+21
19.4
12.88
3.06E+21
17.8
11.83
7.01E+20
10.9
7.24
7.05E+20
10.9
7.25
1.68E+22
6.7
26.24
2.52E+22
8.3
1.88E+21
10.5
12.08
4.04E+20
6.3
7.24
2.20E+21
3.1
31.57
4.48E+21
4.7
5.13E+20
8.5
7.00
2.28E+21
14.0
11.52
5.90E+20
8.9
7.34
5.34E+21
4.4
37.87
8.73E+21
5.9

Ar
2.8
2.6
1.5
1.5
4.0
2.7
1.5
4.6
1.5
2.6
1.5
5.0

大 小
大 小

破壊開始点
アスペリティ（深い位置）
アスペリティ（浅い位置）
都道府県庁所在地

図 3.2.2

図 3.2.3

震源モデル（日本海溝・岩手県沖～青森県沖）（日本海溝・3.11 タイプ）

破壊開始点
アスペリティ（深い位置）
小 アスペリティ（中位）小 アスペリティ（浅い位置）
アスペリティ（内閣府）
都道府県庁所在地

大 小
大大

図 3.2.4

震源モデル

震源モデル（南海トラフ）- 89 -

札幌
苫小牧

大 小
大 小

図 3.2.5

Slip Vel. (m/s)

4

破壊開始点
アスペリティ（深い位置）
アスペリティ（浅い位置）

震源モデル（千島海溝[6]・昨年度計算）

三角形①（面積A1）

3

Ar = A2 / A1
Tr = τ2 / τ1
すべり速度時間関数
（①＋②）

2

三角形②
（面積A2）

1
0

τ1

τ2

0

2

図 3.2.6

すべり速度時間関数（山本・吉村[15]）4Time (s)- 90 -

6

8

10

(a) 3 秒(b) 5 秒(c) 7 秒(d) 10 秒

図 3.2.7

擬似速度応答スペクトル（日本海溝（青森県沖～岩手県沖）・EW）- 91 -

図 3.2.8(a) 3 秒(b) 5 秒(c) 7 秒(d) 10 秒

擬似速度応答スペクトルのばらつき（日本海溝（青森県沖～岩手県沖）・EW）- 92 -

(a) 3 秒(b) 5 秒(c) 7 秒(d) 10 秒

図 3.2.9 擬似速度応答スペクトル（日本海溝（岩手県沖～茨城県沖）・EW）- 93 -

図 3.2.10(a) 3 秒(b) 5 秒(c) 7 秒(d) 10 秒

擬似速度応答スペクトルのばらつき（日本海溝（岩手県沖～茨城県沖）・EW）- 94 -

(a) 3 秒(b) 5 秒(c) 7 秒(d) 10 秒

図 3.2.11

擬似速度応答スペクトル（南海トラフ・EW）(a) 3 秒(b) 5 秒(c) 7 秒(d) 10 秒

図 3.2.12

擬似速度応答スペクトルのばらつき（南海トラフ・EW）- 95 -

(a) 3 秒(b) 5 秒(c) 7 秒(d) 10 秒：3.11 の観測記録の擬似速度応答スペクトル
図 3.2.13

擬似速度応答スペクトル（日本海溝（岩手県沖～茨城県沖）・川辺ほか[13]の 3.11 モデル・EW）- 96 -

参考文献[1] 地震調査研究推進本部：南海トラフの地震活動の長期評価（第二版）について，2014.
[2] 五十嵐さやか・坂本成弘・西田明美・村松健・高田毅士：震源特性の不確定性が地震動強さや建物応答に及ぼす影響，日本建築学会構造系論文集，第 81 巻，第 721 号，pp.425-435，2016.
[3] （一社）日本建築構造技術者協会関西支部 大阪府域内陸直下型地震に対する建築設計用地震動および設計法に関する研究会：大阪府域内陸直下型地震に対する建築設計用地震動および耐震設計指針，
2015.
[4] 地震調査研究推進本部：長周期地震動評価 2016 年試作版－相模トラフ巨大地震の検討－，2016.
[5] 西本昌・山本優・縣亮一郎・藤田航平・市村強・堀高峰：震源特性の不確実性が長周期地震動に与える影響の分析 －2016 年熊本地震前震を対象とした検討－，大成建設技術センター報，第 54 号，31，2021.
[6] 西本昌・谷翼・山本優・縣亮一郎・堀高峰：多様なシナリオを考慮した長周期地震動評価と建物応答解析 千島海溝沿いの巨大地震を対象とした検討，大成建設技術センター報，第 55 号，12，2022.
[7] 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, pp.1526, 2014.
[8] T. Ichimura, M. Hori and J. Bielak：A hybrid multiresolution meshing technique for finite
element three dimensional earthquake ground motion modelling in basins including topography,
Geophysical Journal International, Vol.177，Issue 3, pp.1221-1232, 2009.
[9] 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 to enhance comprehensive earthquake simulation, Proceedingsof the International Conference on High Performance Computing, Networking, Storage andAnalysis, SC’15, 2015.
[10]防 災 科 学 技 術 研 究 所 ： 地 震 ハ ザ ー ド ス テ ー シ ョ ン J-SHIS ， オ ン ラ イ ン （ https://www.jshis.bosai.go.jp/，2022/08/04 参照）[11]地 震 調 査 研 究 推 進 本 部 ： 全 国 1 次 地 下 構 造 モ デ ル （ 暫 定 版 ）， オ ン ラ イ ン（https://www.jishin.go.jp/main/chousa/12_choshuki/choshuki2012_a2.pdf，2022/08/04 参照）[12]地 震 調 査 研 究 推 進 本 部 ： 海 溝 型 地 震 の 長 期 評 価 ， オ ン ラ イ ン（ https://www.jishin.go.jp/evaluation/long_term_evaluation/subduction_fault/ ， 2022/08/04
参照）[13]川辺秀憲・釜江克宏：2011 年東北地方太平洋沖地震の震源のモデル化，日本地震工学会論文集，第13 巻，第 2 号（特集号），2013.
[14]地震調査研究推進本部：震源断層を特定した地震の強震動予測手法（「レシピ」），2020.
[15]山本優・吉村智昭：3 次元大規模 FEM による東海・東南海・南海連動地震の長周期地震動シミュレーション，日本建築学会構造系論文集，第 77 巻，第 677 号，pp.1055-1064，2012.

- 97 -

（４）プロジェクトの総合的推進
プロジェクトの総合的推進では、1) 参画各機関の協力・連携を密にし、円滑にプロジェクトを運営していくため、進捗管理等のための会合等を開催すること、2) 外部有識者を招聘した諮問委員会から、本プロジェクトの進捗状況の把握・評価・改善提言・指導等を受けること、を目的としている。また、プロジェクトで得られた成果については、広報や研究活動を通じて積極的に公表するとともに、成果報告会を開催することとし、「富岳」成果創出加速プログラム事業全体のアウトリーチにも協力する。
1)

進捗管理等のための会合に関しては、昨年度までより一層に効率性を重視した。プロジェクトに参加する研究者全体の連携を密とすることは必須であるが、運営に過度の労力を割くことは決して得策ではないからである。昨年と同じく、進捗管理等のための会合を定期的に開催する方針を維持することとし、全 3 回の全体会議を開催した。
（開催日：令和 4 年 5 月 20 日、9 月 30 日、令和 5
年 1 月 13 日）。昨年度に引き続き、実施機関が開催し、協力機関を交えて進捗状況や今後の予定
等の情報を共有する形で密な連携を図った。今年度は課題代表者を中心に、副代表者や各協力機関のメンバーを交えた小規模な打合せを密にすることで、全体会議の開催回数を少なくし効率化を
図った。

2)

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

東京大学名誉教授

委員
渦岡 良介

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

小野 謙二

九州大学情報基盤研究開発センター 応用データ科学研究部門 教授

山中 浩明

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

矢来 博司

国土地理院 測地観測センター 地震調査官

委員長・委員の専門は、観測地震学・地盤工学・地震工学・数値流体力学・計算科学・物理探査工学・
強震動地震学・固体地球科学である。また、本課題で構築される統合的予測システムを社会実装する際、最大のユーザとなることが想定される国レベルの地震防災・減災の有識者でもある。
第 5 回の諮問委員会は令和 4 年 11 月 24 日に Web 形式で開催された。この委員会では、最終年度計画と進捗状況並びに KPI の進捗確認、加えて今年度が最終年度であることを踏まえ、今後に向けての説明をおこなった。最終年度計画と進捗状況に関しては、本プロジェクトを構成する次の 3 つのサブテーマ
についてそれぞれ計画と進捗状況を説明した。サブテーマ１は、国の被害想定のための地震・津波災害予測でのアプリケーション活用として、国（内閣府防災担当）の長周期地震動評価に関する実施支援を海洋研究開発機構が担当している。サブテーマ２は地震に関する災害被害予測のための大規模アプリケーシ
ョンの開発であり、東京大学地震研究所の市村副代表が中心となって進めている。サブテーマ３は、開発した計算アプリケーション群を実務で活用するための環境整備として、企業がそれぞれの企業ごとに具- 98 -

体的な問題設定で「富岳」ならではの地震動・地盤増幅計算を実施するというものである。各テーマ共順
調に進捗しており、特に大きな成果として、市村副代表等の論文 "Extreme-scale earthquake simulation withuncertainty quantification" が国際会議「SC22」(Supercomputing Conference)にてスーパーコンピュータの世
界で最も権威ある賞のひとつであるゴードン・ベル賞のファイナリストに選出されたことについて紹介
し、委員からは最先端の成果が出ているとの高い評価を受けた。KPI に関しては、次の様な内容の説明をおこなった。1) 国の想定のための地震動計算に関しては、まず、大規模計算の部分において、計算科学としての目標については昨年度の一気通貫のシミュレーションで既に達成済みであり、成果に関する目
標については、国の長周期地震の結果と同等のものを公表し、今年度も内閣府防災担当への協力を続けている。計算に関しては、計算の妥当性確認等のコンサルティング対応に加え、計算環境等の提供を実施した。2) 企業による地震動計算の活用に関しては、連携企業に対し、国（内閣府防災担当）と同等の計算が出来るようなサポートを続け、ユーザビリティの改善をおこなった。関東平野、千島海溝など現実的な問題設定での検討を進めており、今年度達成の見込みである。また今年度は、地震動計算を活用した連携
企業による４件の学会発表があり、今後も学会や論文等の発表につながることが期待される。3) 企業による非線形地盤増幅計算の活用に関しては、「富岳」全系での超大規模計算により、昨年度中に達成済み。
4) 粘弾性応答の参照解作成に関しては、収束性などの品質保証を定量的に行うため、南海トラフを対象にした非常に詳細な構造モデル(0.55 兆自由度)での粘弾性応答計算を東京大学地震研究所にて実施した。
これにより粘弾性の応答の計算は、断層の近傍まで含めて収束した形での計算ができてきていると考えられる。5) 地震波動場計算の参照解作成に関しては、地震波動場計算で品質保証することが目的であり、関東における巨大地震による長周期地震動計算において、2 秒(内閣府での対象周期)までのメッシュサイ
ズを変えた計算による精度確認を今年度中に達成する見込みである。委員からは、「「富岳」を使うことで、これまでには出来なかった自由度の大きな計算が出来るようになったというところは全くその通りなので、もう少し具体的に書いた方がよい。
」「内閣府の評価にどこまで寄与したかを明解にすれば良い。
」
等の助言を受けた。今後に向けて、地震動に関しては、連携企業や内閣府防災担当へのサポートの継続が必要と考えていること、国への取組みとしては地殻変動にも展開していきたいと考えていること、加えて、地震動も含めて大規模計算ができる人材を育てていきたいと考えていることについて言及し、委員から賛同を得られた。人材育成に関しては、固体地球分野でのターゲット人材を明確にすべきとの助言を受けた。
第 6 回の諮問委員会は令和 5 年 2 月 24 日に Web 形式で開催された。この委員会では、最終年度である
今年度の研究成果と KPI 達成状況等を報告し、評価を受けた。KPI の各項目に対する具体的な実施結果と達成状況について表 4.1 に示す。成果報告を受け、委員からは、最先端の計算科学の基礎研究が非常に進展したという事について極めて高い評価を受けた。また、社会実装という観点からは、“学”が開発したコードに対して“産”が機能拡張を連携によって進めていることが非常に頼もしくこの課題の大きな成果だと考える、今後とも引続き“産学官の連携”によってこういった研究が進んでいく事を期待する旨のコメントがあった。結論として、本年度の研究の進捗状況・成果に関して、特段の問題がないことが了解された。

- 99 -

表 4. 1 KPI 達成状況
内容

結果

「富岳」全系で、テラ自由度の 2.5Hz「富岳」全系での超大規模計算により、までの地震動計算を 10Λ4 ステップ行

地盤増幅計算とあわせて令和 3 年度中う。

に達成済み。

国の想定と同等の長周期地震動計算

同等の計算を実施済み。令和 5 年 5 月の結果を公表する。

学会発表をもって達成。

達成

1.国の想定のための地震動計算
計算科学としての目標

成果に関する目標〇〇
2.企業による地震動計算の活用
計算科学としての目標

成果に関する目標「富岳」全系で、テラ自由度の 2.5Hz「富岳」全系での超大規模計算により、までの地震動計算を 10Λ4 ステップ行

地盤増幅計算とあわせて令和 3 年度中う。 (1.と同じ)に達成済み。

国と同等の長周期地震動計算を、建築

関東平野、千島海溝など現実的な問題設系の各企業が独自の問題設定で実施

定で計算を実施、学会でも発表しておする。そのための企業のサポート体制

り、達成したといえる。

〇〇
を海洋研究開発機構に整備する。
3.企業による非線形地盤増幅計算の活用
---

---

独自の問題設定での解析が、東京大学地震研究所との共同研究として進んでおり、さらに独自コードへの大規模ソルバ〇

ーの組み込みまで行っており、想定以上に達成。
4.粘弾性応答の参照解作成
計算科学としての目標

テラ自由度の計算モデルを構築し、数「富岳」全系で曖昧さまで考慮した超大百ステップの粘弾性応答計算を「富規模計算により、令和 4 年度達成。

〇

岳」の全系で実施する。
成果に関する目標

巨大地震の震源域を中心とした 1/8 球

南海トラフを対象にした詳細構造モデの粘弾性応答の参照解とその構造モル(0.55 兆自由度)での粘弾性応答計算

デルを公開する。令和 4 年以降に、連を、令和 4 年度 東京大学地震研究所に携先がこの参照解を活用してグリーて実施したことにより達成。（領域サイン関数ライブラリを構築する。→ 粘

ズ：2496km×2496km×1100km、最小要素

弾性応答計算で品質保証をすること

サイズ：125m、プレート内地震後の 20 年が目的。

間の粘弾性応答）〇5.地震波動場計算の参照解作成
計算科学としての目標

テラ自由度の日本列島モデルを構築
Λ

し、10 4 ステップ以上の弾性波伝播計
算を「富岳」の全系で実施する。

- 100 -

「富岳」全系での超大規模計算により、
令和 3 年度達成済み。

〇

成果に関する目標

巨大地震の震源域を中心とした全国

1.の関東における巨大地震による長周

一次地下構造モデルでの周期 2 秒まで

期地震動計算において、2 秒(内閣府でのの参照解とその構造モデルを公開す

対象周期)までのメッシュサイズを変える。令和 4 年以降に、連携先がこの参た計算による精度確認について令和 4

照解を活用してグリーン関数ライブ

年度にほぼ達成。

〇

ラリを構築する。→地震波動場計算で品質保証をすることが目的
達成欄の○は達成済みであることを示す

プロジェクトで得られた成果については、広報や研究活動を通じて積極的に公表し、今後の展開に資することを継続した。具体的には、国内外の各種学会・研究集会等での発表、および、それら情報を取りまとめた本プロジェクトの和文・英文ホームページ（図 4.1）を介した公表である。ホームページ（和文 https://www.jamstec.go.jp/fugaku-earthq/ja/ 英文 https://www.jamstec.go.jp/fugaku-earthq/en/）では見やすさを考慮し、プロジェクトの最新情報をトップページに時系列で掲載するとともに、イベント、研究成果の別で掲載する形式を維持した。加えて今年度は、イベント、研究成果等のページにて年度毎の表示切替えを可能とすることで閲覧時の視認性に改善を加えた。(図 4.2)
さらに、令和４年度成果報告会(令和 5 年 2 月 24 日 Web 形式)を開催した。開催形式について、昨年
度より引続き Zoom Webinar を利用した。参加者は全体で 106 名だった。成果報告会では本年度の成果に係る報告の後、パネルディスカッションとして連携機関・企業より代表して 6 名にご登壇いただき、本課題代表・副代表を交え、本課題への参画前後の変化や今後への期待について、幅広く闊達な討論が交わされた。
「富岳」成果創出加速プログラム事業全体のアウトリーチへの協力としては、各種学会・研究集会での発表時や、Web マガジン等のメディア取材時に、「富岳」成果創出加速プログラムや、本プロジェクトの紹介、本プロジェクトを通して得られた成果の発表をおこなった。また、イベントを通した協力としては、令和 4 年 5 月に開催された JpGU Meeting 2022 における、堀課題代表を座長とするセッションや、令和 4 年 7 月に開催された NWC CAE World 2022 での同じく堀課題代表による基調講演、令和 4 年 10 月に開催された「第５回 HPCI コンソーシアムシンポジウム」での市村副代表による基調講演、同じく市村副代表による、令和 4 年 12 月に開催された『第２回 スーパーコンピュータ「富岳」成果創出加速プログラ
ムシンポジウム 「富岳百景」』での特別講演、および令和 5 年 1 月に開催された『「富岳」EXPANDS～可能性を拡張する～ 』での講演等が挙げられる。また、本プロジェクトのホームページ上でもイベント告知等の広報活動をおこなった。
計算機環境については、海洋研究開発機構に整備した計算機環境により、重点課題で開発したアプリケーション群や各関係機関で得られた貴重なデータならびに、本年度の業務で得られた成果等の保管・管理がなされている。研究を推進するための計算機環境の維持・整備として、これら計算機の定常的な監視や定期保守を実施した。

- 101 -

図 4. 1 和文・英文ホームページイメージ

図 4. 2 イベント、研究成果ページの年度切替え表示- 102 -

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

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

開催日
令和 4 年 5 月 20 日
令和 4 年 9 月 30 日
令和 5 年 1 月 13 日

委員会

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

令和 4 年 11 月 24 日
令和 5 年 2 月 24 日

成果報告会

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

令和 5 年 2 月 24 日

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

担当機関

担当責任者

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

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

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

堀

アプリケーション活用

海域地震火山部門（２）地震に関する災害被害

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

教授

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

東京大学地震研究所

市村 強ーションの開発

計算地球科学研究センター（３）開発した計算アプリケ

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

センター長ーション群を実務で活用する

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

堀

ための環境整備

海域地震火山部門（４）プロジェクトの総合的神奈川県横浜市金沢区昭和町3173-25

センター長

推進

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

堀

海域地震火山部門- 103 -

高峰

高峰

高峰

ูῧ㸯Ꮫ఍➼Ⓨ⾲ᐇ⦼㸬Ꮫ఍➼࡟࠾ࡅࡿཱྀ㢌࣭࣏ࢫࢱ࣮Ⓨ⾲
Ⓨ⾲ࡋࡓᡂᯝ㸦Ⓨ⾲㢟┠ࠊཱྀ㢌࣭࣏ࢫࢱ࣮Ⓨ⾲ࡢู㸧Ⓨ⾲⪅ẶྡⓎ⾲ࡋࡓሙᡤ
㸦Ꮫ఍➼ྡ㸧

Ⓨ⾲ࡋࡓ᫬ᮇ

ᅜෆ࣭እࡢูᐩᓅࢆ⏝࠸ࡓ኱つᶍᆅ㟈ືィ⟬࡟ྥࡅࡓ9HULILFDWLRQᆅᙧ࡜┅ᆅࢆక࠺ᵓ㐀࡛ࡢ࣋ࣥࢳ࣐࣮ࢡࢸࢫࢺ

ᑠᯘᗈ᫂Ώ㎶ဴྐ➟ᯇ೺ኴ㑻ຍ⸨◊୍⦩ு᪥ᮏᆅ⌫ᝨᫍ⛉Ꮫ㐃ྜᖺ኱఍
୍㑻ᇼ㧗ᓠ



ᅜෆࢫ࣮ࣃ࣮ࢥࣥࣆ࣮ࣗࢱᐩᓅࢆ⏝࠸ࡓ༓ᓥᾏ⁁ἢ࠸ࡢᕧ኱ᆅ㟈ࡢ㛗࿘ᮇᆅ㟈ືࡢホ౯すᮏᫀᒣᮏඃ⦩ு୍㑻ᇼ㧗ᓠ᪥ᮏᆅ⌫ᝨᫍ⛉Ꮫ㐃ྜᖺ኱఍

ᅜෆ

9LVXDOL]DWLRQ0HWKRG8VLQJ3ULQFLSDO$[LV/LQHVDQG,WV$SSOLFDWLRQWR6WUHVV)LHOGGXHWR)DXOW6OLSᘅ㒊⣪ஓᏊࠊ⦩ு୍㑻ࠊᑠᅧ೺஧ࠊᇼ㧗ᓠࠊᇼ᐀ᮁ

᪥ᮏᆅ⌫ᝨᫍ⛉Ꮫ㐃ྜᖺ኱఍

ᐩᓅࢆ⏝࠸ࡓ኱つᶍᆅ㟈ືィ⟬࡟ྥࡅࡓ9HULILFDWLRQ㸫⌧ᐇⓗ࡞ᆅ┙ᵓ㐀࡛ࡢ࣋ࣥࢳ
࣐࣮ࢡࢸࢫࢺ㸫ᑠᯘᗈ᫂Ώ㎶ဴྐ➟ᯇ೺ኴ㑻ຍ⸨◊୍⦩ுᖺᗘ᪥ᮏᘓ⠏Ꮫ఍኱఍୍㑻ᇼ㧗ᓠ

」ᩘࡢⓎ⏕ࢩࢼࣜ࢜ࢆ⪃៖ࡋࡓ㛗࿘ᮇᆅ㟈ືホ౯㸫༓ᓥᾏ⁁ἢ࠸ࡢᕧ኱ᆅ㟈ࢆᑐ㇟࡜ࡋࡓ᳨ウ㸫すᮏᫀᒣᮏඃ⦩ு୍㑻ᇼ㧗ᓠᖺᗘ᪥ᮏᘓ⠏Ꮫ఍኱఍*38ࢆ౑⏝ࡋࡓ㧗㏿࡞ḟඖᆅ㟈᫬ᾮ≧໬ゎᯒ࡜ࡑࡢࢧࣟࢤ࣮ࢺࣔࢹࣝ

᪥ୗ㒊ுኴࠊᕷᮧᙉࠊ⸨⏣⯟ᖹࠊᇼ᐀ᮁࠊ
0DGGHJHGDUD/DOLWKᅵᮌᏛ఍➨ᅇᛂ⏝ຊᏛࢩ࣏ࣥࢪ࣒࢘
ᖺ᭶᪥㹼᪥ࣁ࢖ࣈࣜࢵࢻ㛤ദ࢜ࣥࣛ࢖ࣥཧຍ



ᅜෆ

$SSOLFDWLRQRI0HWD0RGHOLQJ7KHRU\WR7KLQ&XUYHG%HDP8VLQJ&XUYLOLQHDU&RRUGLQDWH6\VWHPDQG3HUWXUEDWLRQ([SDQVLRQ

0XQHR+RUL.RKHL)XMLWD

ᅵᮌᏛ఍➨ᅇᛂ⏝ຊᏛࢩ࣏ࣥࢪ࣒࢘
ᖺ᭶᪥㹼᪥ࣁ࢖ࣈࣜࢵࢻ㛤ദ࢜ࣥࣛ࢖ࣥཧຍ

ᅜෆ㒔ᕷࡢヲ⣽ḟඖᆅ㟈ࢩ࣑࣮ࣗࣞࢩࣙࣥࡢࡓࡵࡢ㧗㏿࣭ࢫࢣ࣮ࣛࣈࣝ࡞᭷㝈せ⣲ἲࡢ㛤Ⓨ⸨⏣⯟ᖹࠊᕷᮧᙉࠊᇼ᐀ᮁࠊ0DGGHJHGDUD/DOLWK➨ᅇィ⟬ᕤᏛㅮ₇఍ᖺ᭶᪥㹼᪥㸦⛅⏣┴⛅⏣ᕷ㸧ࣁ࢖ࣈࣜࢵࢻ㛤ദ࢜ࣥࣛ࢖ࣥཧຍ



ᅜෆ࣓ࢽ̿ࢥ࢔&38ཬࡧ*38ࢆ౑⏝ࡋࡓ᭷㝈せ⣲ἲ࡟ࡼࡿ኱つᶍḟඖᆅ┙ᾮ≧໬ゎᯒ࡜ࡑࡢࢧ ᪥ୗ㒊ுኴࠊᕷᮧᙉࠊ⸨⏣⯟ᖹࠊᇼ᐀ᮁࠊ
ࣟࢤ࣮ࢺࣔࢹࣝ
:LMHUDWKQH/DOLWK➨ᅇィ⟬ᕤᏛㅮ₇఍ᖺ᭶᪥㹼᪥㸦⛅⏣┴⛅⏣ᕷ㸧ࣁ࢖ࣈࣜࢵࢻ㛤ദ࢜ࣥࣛ࢖ࣥཧຍ



ᅜෆ

7KH,QWHUQDWLRQDO&RQIHUHQFHRQ&RPSXWDWLRQDO6FLHQFH,&&6WK,QWHUQDWLRQDO:RUNVKRSRQ$GYDQFHVLQ+LJK3HUIRUPDQFH&RPSXWDWLRQDO(DUWK6FLHQFHV$SSOLFDWLRQV )UDPHZRUNV

ᅜእ

ᅜእ

6&:$&&3'1LQWK:RUNVKRSRQ$FFHOHUDWRU3URJUDPPLQJ
8VLQJ'LUHFWLYHV :$&&3' 

ᅜእ&DOFXODWLRQRI&URVVFRUUHODWLRQ)XQFWLRQ$FFHOHUDWHGE\7HQVRU&RUHVZLWK
7HQVRU)ORDWSUHFLVLRQRQ$PSHUH*38<XPD.LNXFKL.RKHL)XMLWD7VX\RVKL,FKLPXUD0XQHR+RULDQG/DOLWK0DGGHJHGDUD([WUHPH6FDOH(DUWKTXDNH6LPXODWLRQZLWK8QFHUWDLQW\4XDQWLILFDWLRQ㸪6&㸦*RUGRQ%HOO3UL]H)LQDOLVW㸧

7VX\RVKL,FKLPXUD.RKHL)XMLWD5\RWD.XVDNDEH.HQWDUR.R\DPD6RWD0XUDNDPL
<XPD.LNXFKL7DNDQH+RUL0XQHR+RUL

*38$FFHOHUDWHG6SDUVH0DWUL[9HFWRU3URGXFW%DVHGRQ(OHPHQWE\(OHPHQW0HWKRG 5\RWD.XVDNDEH.RKHL)XMLWD7VX\RVKL
IRU8QVWUXFWXUHG)(08VLQJ2SHQ$&&,FKLPXUD0XQHR+RUL0DGGHJHGHUD/DOLWK

6&,QWHUQDWLRQDO&RQIHUHQFHIRU+LJK3HUIRUPDQFH&RPSXWLQJ1HWZRUNLQJ6WRUDJHDQG$QDO\VLV*RUGRQ%HOO3UL]H)LQDOLVW





ᅜෆ

ᅜෆ

ᅜෆ

6FDODEOH)LQLWH(OHPHQW9LVFRHODVWLF&UXVWDO'HIRUPDWLRQ$QDO\VLV$FFHOHUDWHG
ZLWK'DWD'ULYHQ0HWKRG.RKHL)XMLWD6RWD0XUDNDPL7VX\RVKL,FKLPXUD7DNDQH+RUL0XQHR+RUL
0DGGHJHGDUD/DOLWKDQG1DRQRUL8HGD

6&6FDO$+WK:RUNVKRSRQ/DWHVW$GYDQFHVLQ
6FDODEOH$OJRULWKPVIRU/DUJH6FDOH+HWHURJHQHRXV
6\VWHPV



ᅜእᵓ㐀ࡢศ࠿ࡽ࡞࠸ᆅẆࡢኚᙧࢆண

ᕷᮧࠉᙉ

5,67+3&,࣐࢞ࢪࣥࠕᐩᓅⓒᬒࠖYRO

ᅜෆ༡ᾏࢺࣛࣇࡢᕧ኱ᆅ㟈࡟ࡼࡿ⢓ᙎᛶኚᙧ࡟ᑐࡍࡿࣜࢯࢫࣇ࢙࢔࢔ࢭࣀࢫࣇ࢙࢔ቃ⏺ࡢప ᮧୖ㣁ኴࠊᶫ㛫᫛ᚨࠊ㣤἟༟ྐࠊ⸨⏣⯟ᖹࠊᕷ᪥ᮏᆅ㟈Ꮫ఍ᖺᗘ⛅Ꮨ኱఍ࠊ
⢓ᛶ㡿ᇦࡢᙳ㡪ᮧᙉࠊᇼ㧗ᓠ

ᅜෆ㧗ヲ⣽୕ḟඖᆅ┙㟈ືゎᯒᡭἲࡢ㛤Ⓨ࡜ࠊᆅ┙ṍࡢ࣮࣒࣎ࣜࣗࢹ࣮ࢱྍど໬࡟㛵ࡍࡿᇶ
♏᳨ウ

⸨⏣⯟ᖹ

᪥ᮏᆅ㟈ᕤᏛ఍ࠉ-$(($QQXDO0HHWLQJᖺ᭶᪥㹼᭶᪥ࠊᮐᖠ࣭໭ᾏ㐨኱Ꮫ

ᅜෆ≉ูㅮ₇ࠉࠕᐩᓅࠖ࡟ࡼࡿᆅ㟈ࢩ࣑࣮ࣗࣞࢩࣙࣥ

ᕷᮧࠉᙉᩥ㒊⛉Ꮫ┬ࢫ࣮ࣃ࣮ࢥࣥࣆ࣮ࣗࢱࠕᐩᓅࠖᡂᯝ๰ฟຍ㏿ࣉࣟࢢ࣒ࣛࠉ➨ᅇᡂᯝ๰ฟຍ㏿ࣉࣟࢢ࣒ࣛࢩ࣏ࣥࢪ࣒࢘ࠕᐩᓅ
ⓒᬒࠖ㺎ࠕᐩᓅࠖᖺࡢ㌶㊧ࠊᮾி㒔ࠊ㸦࢜ࣥࣛ࢖ࣥ㛤ദ㸧

ᅜෆࢩ࣑࣮ࣗࣞࢩࣙࣥ࡜ࢹ࣮ࢱࡢ⼥ྜ࡟ࡘ࠸࡚ᆅ㟈ࢆ౛࡟ᕷᮧࠉᙉࢩ࣏ࣥࢪ࣒࢘ࠕᐩᓅࠖ(;3$1'6㹼ྍ⬟ᛶࢆᣑᙇࡍࡿ㹼ࣃࢿࣝࢹ࢕ࢫ࢝ࢵࢩࣙࣥࠕᐩᓅࠖ12:㸸ࠕᐩᓅࠖࢆ⏝࠸ࡓ◊✲ᡂᯝࠊᮾி㒔ࠊ㸦ࣁ࢖ࣈࣜࢵࢻ㛤ദ㸧

ᅜෆ

7KHWK5&&6,QWHUQDWLRQDO6\PSRVLXP㸦ࣁ࢖ࣈࣜࢵࢻ㛤ദ㸧

ᅜෆࡍࡿᖖ㆑ࢆ㉸࠼ࡓ኱つᶍᆅ㟈ࢩ࣑࣮ࣗࣞࢩࣙࣥ(QKDQFHG(DUWKTXDNH6LPXODWLRQZLWK'DWD'ULYHQ0HWKRGVDQG6WRFKDVWLF$QDO\VLV .RKHL)XMLWD- 104 -

㸰㸬Ꮫ఍ㄅ࣭㞧ㄅ➼࡟࠾ࡅࡿㄽᩥᥖ㍕Ⓨ⾲ࡋࡓሙᡤ
㸦Ꮫ఍ㄅ࣭㞧ㄅ➼ྡ㸧Ⓨ⾲ࡋࡓ᫬ᮇ

ᅜෆ࣭እࡢูከᵝ࡞ࢩࢼࣜ࢜ࢆ⪃៖ࡋࡓ㛗࿘ᮇᆅ㟈ືホ౯࡜ᘓ≀ᛂ⟅ゎᯒ༓ᓥᾏ⁁ἢ࠸ࡢᕧ኱ᆅ㟈ࢆ すᮏᫀࠊ㇂⩼ᒣᮏඃࠊ⦩ு୍㑻ࠊᇼ㧗ᑐ㇟࡜ࡋࡓ᳨ウᓠ኱ᡂᘓタᢏ⾡ࢭࣥࢱ࣮ሗ12
KWWSVZZZWDLVHLFRMSWDFWWUSDSHU$BSGIᅜෆ

&RUUHFWLRQWR+LJKILGHOLW\HODVWLF*UHHQ̓VIXQFWLRQVIRUVXEGXFWLRQ]RQH
PRGHOVFRQVLVWHQWZLWKWKHJOREDOVWDQGDUGJHRGHWLFUHIHUHQFHV\VWHP

7DNDQH+RUL5\RLFKLUR$JDWD7VX\RVKL,FKLPXUD.RKHL)XMLWD7DNXPD<DPDJXFKL7DNHVKL,LQXPD($57+3/$1(76$1'63$&($UWLFOHQXPEHU '2,KWWSVGRLRUJV[

ᅜእ,PSDFWRI$PELJXLW\RI3K\VLFDO3URSHUWLHVRI7KUHH'LPHQVLRQDO&UXVWDO
6WUXFWXUH0RGHORQ&RVHLVPLF6OLSDQG,QWHUVHLVPLF6OLS'HILFLWLQWKH1DQNDL
7URXJK5HJLRQ

6RWD0XUDNDPL7VX\RVKL,FKLPXUD.RKHL)XMLWD7DNDQH+RUL<XVDNX2KWD*HR+D]DUGV$GYDQFHG1XPHULFDO6LPXODWLRQIRU(DUWKTXDNH+D]DUGVDQG'LVDVWHUV

ᅜእ$1RQ6LJQDOL]HG-XQFWLRQ0RGHOIRU$JHQW%DVHG6LPXODWLRQVRI&DU̽3HGHVWULDQ
0RGH0DVV(YDFXDWLRQV

0DGGHJHGDUD/DOLWK:DVXZDW3HWSUDNRE
0XQHR+RUL7VX\RVKL,FKLPXUD.RKHL)XMLWD*HR+D]DUGV$GYDQFHG1XPHULFDO6LPXODWLRQIRU(DUWKTXDNH+D]DUGVDQG'LVDVWHUV

ᅜእ

6FDODEOHODUJHVFDOHPXOWLSK\VLFVHDUWKTXDNHVLPXODWLRQRQPXOWLSOH*38VZLWK
VWDELOL]DWLRQ

5\RWD.XVDNDEH7VX\RVKL,FKLPXUD.RKHL)XMLWD0XQHR+RUL
/DOLWK:LMHUDWKQH-RXUQDORI&RPSXWDWLRQDO6FLHQFH
YRO

ᅜእ&DOFXODWLRQRI&URVVFRUUHODWLRQ)XQFWLRQ$FFHOHUDWHGE\7HQVRU&RUHVZLWK
7HQVRU)ORDWSUHFLVLRQRQ$PSHUH*38,&&6

<XPD.LNXFKL.RKHL)XMLWD7VX\RVKL,FKLPXUD0XQHR+RUL0DGGHJHGDUD/DOLWK

&RPSXWDWLRQDO6FLHQFH,&&6



ᅜእィ⟬ᕤᏛㅮ₇఍ㄽᩥ㞟3URFHHGLQJVRIWKH&RQIHUHQFHRQ&RPSXWDWLRQDO(QJLQHHULQJDQG6FLHQFH᪥ᮏィ⟬ᕤᏛ఍⦅



ᅜෆ

ᅜእ

ᥖ㍕ࡋࡓㄽᩥ㸦Ⓨ⾲㢟┠㸧Ⓨ⾲⪅Ặྡ࣓ࢽ࣮ࢥ࢔&38ཬࡧ*38ࢆ౑⏝ࡋࡓ᭷㝈せ⣲ἲ࡟ࡼࡿ኱つᶍḟඖᆅ┙ᾮ≧໬ゎᯒ࡜ࡑࡢࢧ
᪥ୗ㒊ுኴᕷᮧᙉ⸨⏣⯟ᖹᇼ᐀ᮁ
ࣟࢤ࣮ࢺࣔࢹࣝ
)DVWDQG/DUJHVFDOH7KUHHGLPHQVLRQDO6RLO/LTXHIDFWLRQ6LPXODWLRQZLWK)LQLWH /DOLWK:LMHUDWKQH(OHPHQW0HWKRG8VLQJ0DQ\FRUH&38VDQG*38VDQG,WV6XUURJDWH0RGHO
7VX\RVKL,FKLPXUD.RKHL)XMLWD5\RWD([WUHPH6FDOH(DUWKTXDNH6LPXODWLRQZLWK8QFHUWDLQW\4XDQWLILFDWLRQ㸪6&
.XVDNDEH.HQWDUR.R\DPD6RWD
㸦*RUGRQ%HOO3UL]H)LQDOLVW㸧0XUDNDPL<XPD.LNXFKL7DNDQH+RUL0XQHR+RUL+LNDUX,QRXH7DNDIXPL1RVH

6&,QWHUQDWLRQDO&RQIHUHQFHIRU+LJK3HUIRUPDQFH&RPSXWLQJ1HWZRUNLQJ6WRUDJHDQG$QDO\VLV*RUGRQ%HOO3UL]H)LQDOLVW

*38$FFHOHUDWHG6SDUVH0DWUL[9HFWRU3URGXFW%DVHGRQ(OHPHQWE\(OHPHQW0HWKRG 5\RWD.XVDNDEH.RKHL)XMLWD7VX\RVKL
IRU8QVWUXFWXUHG)(08VLQJ2SHQ$&&,FKLPXUD0XQHR+RUL0DGGHJHGHUD/DOLWK

6&:$&&3'1LQWK:RUNVKRSRQ$FFHOHUDWRU3URJUDPPLQJ
8VLQJ'LUHFWLYHV :$&&3' 



ᅜእ

6&6FDO$+WK:RUNVKRSRQ/DWHVW$GYDQFHVLQ
6FDODEOH$OJRULWKPVIRU/DUJH6FDOH+HWHURJHQHRXV
6\VWHPV



ᅜእ

6FDODEOH)LQLWH(OHPHQW9LVFRHODVWLF&UXVWDO'HIRUPDWLRQ$QDO\VLV$FFHOHUDWHG
ZLWK'DWD'ULYHQ0HWKRG.RKHL)XMLWD6RWD0XUDNDPL7VX\RVKL,FKLPXUD7DNDQH+RUL0XQHR+RUL/DOLWK0DGGHJHGDUDDQG1DRQRUL8HGD- 105 -