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

成果報告書

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

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

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

学会等発表実績

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

１．補助事業の目的
本研究では、国難とされる首都直下地震や南海トラフ地震に備えるために必要となる、地震災害の定量的評価手法を、国が実施する地震動・地盤増幅評価に実装することを目指す。そのために、「富岳」での超大規模計算を実施するための統合的予測システムを構築する。それとともに、地震発生予測に工学的な品質保証の概念を導入して政策判断等に耐える基礎を築く。そのうえで、従来、重要構造物の設計に際し、同じ想定地震に対して個別に地震動・地盤増幅の評価を行ってきた建築・土木系の企業が、国の被害想定のために実装・使用されるものと同等の計算ツールやデータを実務で活用するための仕組みを構築する。
２．令和２年度（報告年度）の実施内容
２－１．当該年度（令和２年度）の事業実施計画（１）国の被害想定のための地震・津波災害予測でのアプリケーション活用
重点課題から開発してきた大規模有限要素コードを用いた長周期地震動の計算結果が、国の被害想定に向けた委員会でこれまで検討対象となってきた。本年度は、同様な計算での精度検証作業やそのために必要な計算アプリケーションの修正を行う。
「富岳」を用いた計算が予定される非線形地盤増幅評価に関して
は、大規模有限要素コードによる計算を国と連携して行うための準備を進める。また、津波評価に必要な、地震による海底地殻変動計算に関しても、従来手法から有限要素法へ転換して計算した結果が国の被害想定に向けた委員会で津波評価の一部分としてこれまで検討対象となってきた。本年度は、大規模有限要素コードによるさらに詳細な計算を連携して進める。
（２）地震に関する災害被害予測のための大規模アプリケーションの開発「富岳」の性能を引き出すように、計算科学・計算機科学の最先端技術を駆使して、地震に関する災害
被害予測のための大規模アプリケーションを改良する。まず、地震に関する災害被害予測のためのアプリ
ケーション群の基幹数値解析コンポーネントである、地震・地震動・地殻変動計算等を行う大規模数値解
析用コード群に対し、ソルバーの前処理等の改良を行い、「富岳」実機上でより良好な性能を発揮するように開発を行う。また、大規模数値解析コード群により解析する、全球モデルを視野に入れた超大規模解析モデルの構築手法を解析手法の特性を考慮しつつ開発する。本年度は、より効率的な大規模数値解析を可能とするよう、構造モデルの詳細度にあわせて超大規模モデルのデータ構造を変更した手法の開発を行う。また、構築手法の妥当性を検討するとともに、構築したモデルを用いて大規模数値解析コード群により地震・地震動・地殻変動等の数値解析を行う。
（３）開発した計算アプリケーション群を実務で活用するための環境整備
国の被害想定算出に使われるものと同等の長周期地震動計算・非線形地盤増幅計算を、建築・土木系の各企業が独自の問題設定で実施するためのサポート体制を、海洋研究開発機構が中心となって整備する。
本年度は、プロジェクトに連携機関として入っていただいている企業の方々が計算アプリケーション群を使用するための環境構築を進める。

1

（４）プロジェクトの総合的推進
プロジェクトの総合的推進にあたっては、参画各機関の協力・連携を密にし、円滑にプロジェクトを運営していくため、進捗管理等のための会合等を開催するとともに、外部有識者を招聘した諮問委員会から、本プロジェクトの進捗状況の把握・評価・改善提言・指導等を受ける。また、プロジェクトで得られた成果については、広報や研究活動を通じて積極的に公表する。さらに、重点課題で開発したアプリケー
ション群や各関係機関で得られた貴重なデータならびに、本年度の業務で得られる成果等を保管・管理
し、研究を推進するための計算機環境を維持・整備する。
２－２．実施内容（成果）（１）国の被害想定のための地震・津波災害予測でのアプリケーション活用
課題参加者らは、先行課題であるポスト「京」重点課題から継続して、「富岳」の計算性能を最大限活用できるような有限要素計算アプリケーションの開発に取り組んできた。重点課題において開発された地震
動解析プログラム「E-wave FEM」
、非線形地盤増幅計算プログラム「STRIKE」、地殻変動計算プログラ
ム「E-cycle FEM」は、国の被害想定における計算手法高度化のニーズと合致しており、国の被害想定への適用を目指した検討にすでに組み込まれている。これらに基づき、「富岳」でしかできない規模の地震
動・地盤増幅・地殻変動計算による被害予測システムを構築し、国の被害想定の過程に実装することが本
課題における大きな目標となる。本項目では、まず①において、E-wave FEM を用いた長周期地震動の被害想定に向けた検討について説明する。ここでは、E-wave FEM が国の被害想定のための長周期地震動評
価に使用されるために必要な、計算の精度検証や必要な機能追加を行ったことと、国の長周期地震動評価のための計算を富岳で実施するための準備が整ったことについて述べる。②においては、E-cycle FEM を用いた津波評価のベースとなりえる地殻変動計算に関する検討について述べる。
①

E-wave FEM の精度検証と機能追加将来の発生が懸念される海溝型巨大地震による長周期地震動の推定は、都市部の地震被害対策を進めていくうえで重要である。そのために、適切に設定された震源断層モデルに対し複雑な地盤構造における地震波
動場を数値計算することが必要である。国の被害想定では従来、有限差分法による地震動解析手法が使われてきた。一方、有限差分法は構造格子に基づいて計算を行うため、速度コントラストのある層境界面の形状を滑らかに再現できず、結果的に計算精度が悪くなる場合があることが指摘されている[1]。このような場合、有限要素法のような非構造メッシュに基づいた計算手法を用いることが望ましい。有限要素法導入の障壁となっていたのは、計算コストの大きさである。一方で、本課題参加者のうち、東京大学地震研究所のグループが開発してきた地震動の有限要素計算用プログラム E-wave FEM(例えば[2][3])は、四面体二次要素による詳細な形状表現、陰的時間積分による安定な計算、高性能前処理付き反復解法による高速かつスケーラブ
ルな一次方程式求解により、複雑地盤における大規模地震動計算を精度よく、安定に、かつ高速に行うことを可能としている。また、成層構造地盤に最適化した全自動有限要素メッシュ生成ツールも同梱しており、計算モデルの生成に関わる計算コストの問題も解決している。E-wave FEM の高性能なソルバーはこれまで国内外のフラグシップレベルのスパコンにおいて良好な性能を出しており、高性能計算の分野で高い評価を得ている。そのため、国が行うような被害想定においても、有限差分法による地震動解析を E-wave FEM
によるものに置き換え、幅広い問題設定に対する検討をより精度よく行えるようになることが期待されている。本課題においては、E-wave FEM を社会実装により適したものに改善し、国の被害想定のための長周期地震動評価に使用できると判断されるような準備を行うことが重要な目標の一つとされている。そのために、様々な問題設定における精度検証と必要に応じた機能追加を検討することが重要となる。今年度は、解析解が存在する問題における解析解との比較、国の被害想定で用いられた実績のある既往手法との比較などを通2

じ、精度検証による必要機能の洗い出しと、機能追加のための開発を実施した。そのうえで、国の被害想定のための長周期地震動評価に必要なものと同等の地震動計算を行い、E-wave FEM を国の被害想定に適用するための準備ができていることを確認した。なお、本項目に関する数値計算には、特に断りがない場合、東京大学と筑波大学が共同運営する Oakforest-PACS を用いた。
①-(a). E-wave FEM の概要
E-wave FEM で用いられる計算手法について簡潔に述べる。E-wave FEM では弾性波動の支配方程式を有限要素定式化により離散化したシステムを扱う。すなわち(1-1)

𝐊𝐊𝐊𝐊 + 𝐂𝐂𝐂𝐂 + 𝐌𝐌𝐌𝐌 = 𝐟𝐟である。ここで𝐊𝐊, 𝐂𝐂, 𝐌𝐌, 𝐮𝐮, 𝐯𝐯, 𝐚𝐚, 𝐟𝐟はそれぞれ、剛性マトリクス、ダンピングマトリクス、質量マトリクス、変位
ベクトル、速度ベクトル、加速度ベクトル、外力ベクトルである。ダンピングマトリクスはレイリー減衰モデル[4]により設定される。時間ステップ幅に対する無条件安定性を重視し、陰的な時間積分スキームである
Newmark-𝛽𝛽法(𝛽𝛽=1/4, 𝛿𝛿=1/2)を適用することで、上記の方程式を𝐀𝐀𝛿𝛿𝐮𝐮𝑛𝑛 = 𝐛𝐛𝑛𝑛

(1-2)

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

(1-3)

ただし、

𝐀𝐀 =

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

𝐪𝐪𝑛𝑛 = 𝐪𝐪𝑛𝑛−1 + 𝐊𝐊𝛿𝛿𝐮𝐮𝑛𝑛 ,

(1-5)

2
𝛿𝛿𝐮𝐮 ,
Δ𝑡𝑡 𝑛𝑛

(1-7)

(1-6)

𝐮𝐮𝑛𝑛 = 𝐮𝐮𝑛𝑛−1 + 𝛿𝛿𝐮𝐮𝑛𝑛 ,𝐯𝐯𝑛𝑛 = 𝐯𝐯𝑛𝑛−1 +𝐚𝐚𝑛𝑛 = 𝐚𝐚𝑛𝑛−1 −

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

(1-8)

ここで増分形式をとっているのは非線形問題などへの拡張性を考慮したためであり、弾性波動方程式を解く
場合には必ずしも増分形式とする必要はない。毎時間ステップにおいて、巨大な疎行列𝐀𝐀（最も小規模なケースでも 107 オーダーの次元数）を含む線形方程式の求解が必要となる。共役勾配法・マルチグリッド法・精度混合演算を組み合わせた並列ソルバーを用い、線形方程式の求解を行う[3]ことにより、「富岳」や「京」コンピュータなどのスーパーコンピュータを使用して多自由度（107-109 のオーダー、あるいはしばしばそれ以上）の問題を高速に解くことが可能となる。
3

地震動の数値計算は、弾性波動方程式を離散化して解くことに加え、複雑な地盤構造に対する計算モデルの生成、減衰のモデリング、モデル端部における反射波の処理など、様々な要素技術により実現されるものである。以降の節では、
・E-wave FEM の搭載するこれらの要素が国の被害想定のための長周期地震動評価に対して十分に機能しているかを、解析解やこれまで国の被害想定に用いられてきたものと同等の従来手法との比較による精度検証を通じて確認
・不十分な部分について機能追加したことについて説明する。精度検証の対象となる項目として、A.「波動方程式の求解」、B.「計算モデルの生成」、C.「震源断層のモデリング」、D.「モデル端部における反射波の処理」、E.「減衰のモデリング」の 5つを設定する。表(1).1 に、5 項目に対応する、もともとの E-wave FEM に搭載された手法、精度検証に基づいた対応、精度検証・対応についての説明がある節を示す。
表(1).1 E-wave FEM の精度検証と機能追加の概要
精度検証の対象となる

E-wave FEM にもとも精度検証に基（1）①において精度検証・対応事項と搭載された手法

づいた対応

についての説明がある節

A.波動方程式の求解有限要素法（非構造四面

機能追加の必①-(b). 2 層問題における解析解・体二次要素）
・Newmar

要なし

差分法との比較Squash factor①-(b). 2 層問題における解析解・

k-β 法
B.計算モデルの生成

Ichimura et al.(2009)[2]による非構造メッシ と 要 素 再 分 割差分法との比較ュ生成①-(e). 2013 年 4 月 13 日に淡路

手法の導入島で発生した地震に対する地震
動計算
C.震源断層のモデリン

点震源

点震源と節点①-(c). 面的な広がりを持った震グ

節点分割法[5]分割法の結果源モデルに対する地震動計算のには差がある

精度検証

ことに留意
D.モデル端部における粘性境界条件[6]反射波の処理吸 収 帯 [7] の 導①-(d). 表面波が長時間持続す入る水平 2 層問題における解析解との比較

E.減衰のモデリング

レイリー減衰[4]

機能追加の必①-(d). 表面波が長時間持続す

要なし

る水平 2 層問題における解析解との比較①-(b). 2 層問題における解析解・差分法との比較
まず、堆積層と母岩の速度コントラストの大きい堆積盆地での地震動解析を、E-wave FEM を用いて行う。ここで扱う問題は、地表面で得られる波形が主に実体波によって構成されるような問題であり、A.「波動方程式の求解」と B.「計算モデルの生成」についての検証に用いることができる。具体的には、水平 2 層問題における地震波動場解析解との比較による数値検証、単純だが波動場における 3 次元的効果の出やすい
3 次元不均質構造における地震動計算を行う。国の被害想定に使用されるものと同等の有限差分法でも同じ問題を解くことで、E-wave FEM が従来法と同等、もしくはそれ以上の精度でこれらの問題に対する数値解
4を得られることを示す。
比較対象の差分法について簡単に述べる。用いるのは[8]をベースとした計算プログラムであり、国の被害想定における長周期地震動の計算のために十分な機能を備えているとみなされているものである。空間方向には四次近似、時間積分は二次の陽的スキームを用いている。メモリバリアブル法による減衰のモデリングを行っているため、広く用いられる Q 値が周波数に比例する手法とは異なり、周波数に依存しない減衰の導入が可能となっている。長周期地震動計算用の地盤モデルにおいては、周波数によらない一定の減衰率が想定されていることが一般的であるため、そのような減衰率を計算に取り込める点は、本計算プログラムの長所と言える。なお、これは E-wave FEM で採用されるレイリー減衰モデルにも（近似的ではあるが）共通する特徴である。反射波の処理のために吸収境界が導入されている。
扱う複数の問題設定に共通して、E-wave FEM においては、対象とする最大周波数について、1 波長当たり最低 5 要素となるよう計算メッシュを生成する。二次要素を用いているため、1 波長当たり 10 節点に相当する。比較する差分法においては、1 波長当たり最低 10 格子となる格子間隔で計算を行う。解析解との比較を目的とした水平 2 層モデルを用いた問題設定における震源パラメータと物性値を表(1).2 に示す。ここで、堆積層の厚さを 1km とする。1 次元構造問題での解析解の計算には[9]の手法を用いる。対象とする
最大周波数を 1Hz として地震動計算を行う。堆積層の S 波速度が 1000m/s であるため、1Hz の波 1 波長分
当たり 5 要素のメッシュを生成するため、最小要素サイズは 200m となった。図(1).1 に計算波形の一例として、原点における速度波形を示す。ここに示されるように、水平 2 層問題において、解析解、E-wave FEM、差分法の波形は見た目上よく一致する。ここで、波形の一致度を簡便に評価するため、波形の相互相関係数をもとにしたスコア𝐶𝐶 ∗を考える。すなわち、二つの時刻歴波形𝑎𝑎1 (𝑡𝑡), 𝑎𝑎2 (𝑡𝑡)を考えた時に(1-9)

𝐶𝐶 ∗ = 10max[𝐶𝐶(𝑎𝑎1 (𝑡𝑡), 𝑎𝑎2 (𝑡𝑡)),0]

ここで𝐶𝐶(𝑎𝑎1 , 𝑎𝑎2 ) =∫ 𝑎𝑎1 (𝑡𝑡)𝑎𝑎2 (𝑡𝑡)𝑑𝑑𝑑𝑑
2[∫ 𝑎𝑎1 (𝑡𝑡)𝑑𝑑𝑑𝑑]1/2 [∫ 𝑎𝑎22 (𝑡𝑡)𝑑𝑑𝑑𝑑]1/2(1-10)である。𝐶𝐶 ∗ = 10であれば 2 つの波形は同一のものということになる。これは[10]の提案する波形のGOODNESS-OF-FIT を表す 10 のスコアのうち、
特に時刻歴波形の一致度に直接的にかかわるものとして、用いることとした。𝐶𝐶 ∗の計算時に 10 をかけているのは、GOODNESS-OF-FIT が波形の一致度に対して 10

点満点の評価を与えることを目的として設定されていることに由来する（観測波形と計算波形の比較を想定したもので、4 未満は「poor」、4-6 未満は「fair」、6-8 未満は「good」、8 以上は「excellent」となる）。この指標をもとにすると、原点での x 成分と z 成分の計算速度波形について、E-wave FEM と解析解の𝐶𝐶 ∗ス

コアはそれぞれ 9.98 と 9，97、差分法と解析解の𝐶𝐶 ∗スコアはそれぞれ 9.88 と 9.95 となり、両者ともスコアが高いうえ、E-wave FEM と解析解が特によく一致していることがわかる（y 成分を議論していないのは振幅が小さいため）。
次に単純だが波動場における 3 次元的効果の出やすい問題設定として、堆積層と母岩の速度境界が図(1).2のような形状をした問題を設定する。これは「Moon Basin 問題」と呼ばれる設定[11]であり、母岩にあたる

硬質な地盤層からの地震動に対し、堆積層にあたる軟弱地盤において反射波と透過波の混在による複雑な波形が形成されることが知られている。Moon Basin 問題において E-wave FEM と有限差分法の計算波形を比較する。Moon Basin 問題における震源パラメータを図(1).2 に示す。物性値は水平 2 層問題と同様のものを用いる。この地盤構造に対して E-wave FEM と差分法において生成した計算モデルの一部を図(1).3 に示
5す。非構造四面体二次要素を用いる前者と構造格子を用いる後者での、計算モデルの違いが示されている。
具体的には、E-wave FEM において青線で表される真の境界形状を非常によく再現した計算メッシュが生成されているのに対し、差分法においては構造格子に従い境界形状が階段状にモデル化されている。このよ
うな速度境界の詳細な形状再現性は、公表されている地盤構造モデルにおける地震動の数値解を精度よく計算するという verification の観点から、有限要素法を長周期地震動計算に採用するうえでの一つの利点と言える。図(1).4 に、図(1).2 で示される観測点上での計算速度波形を示す。盆地の内部でも両手法の波形が比較的よく一致する地点が多い。一方、盆地の中央部ではややずれが大きく見える点も存在する。図(1).5 に、盆地の中央部に位置する観測点 B・C・D での速度波形と疑似速度応答スペクトル(h=0.05)を示す。観測点
D の y 成分では、E-wave FEM と差分法の速度波形の𝐶𝐶 ∗スコアは 9.27 と計算波形同士の比較としてはやや低くなっており、疑似応答スペクトルで見ても広い周期帯にわたって相対差が 50%程度となっている。Ewave FEM の解のメッシュサイズに対する収束性を確認するため、すでに得られている結果と 1 波長当たり 10 要素という計算条件で計算した結果の比較の例（観測点 D）を図(1).6 に示す。両者は見た目にもよく一致しており、y 成分についての 1 波長当たり 5 要素の計算波形との𝐶𝐶 ∗スコアは 9.97 となった。このように E-wave FEM の計算波形は、計算メッシュを細かくしても数値解がほとんど変わらないため、空間的な離散化が十分であることが示される。このことから、E-wave FEM と差分法の結果のずれの原因は、Moon

Basin の速度境界面の構造格子による離散化の精度にあるものと推測する。この結果は、複雑な形状を効率よく表現できる有限要素法の利点が示された例と考えられる。これより、E-wave FEM が従来法と同等、もしくはそれ以上の精度で対象問題に対する数値解を得られるレベルにあることが示された。
以上の検討の結果をまとめる。水平成層問題では、E-wave FEM、差分法、解析解はよく一致することが確認された。特に E-wave FEM と解析解はよく一致していた。また、Moon Basin 問題で、E-wave FEM と差分法の堆積盆地内での計算速度波形を比較すると、多くの点ではよく一致した。一方で、いくつかの地点
では、E-wave FEM と差分法の時刻歴波形・疑似応答スペクトルの差が少し大きかった。E-wave FEM では、一般的な離散化条件で、メッシュサイズに対してよく収束した数値解が求められていると考えられる。
これより、A.「波動方程式の求解」
・B.「計算モデルの生成」について、E-wave FEM が従来の国の長周期地震動評価に用いられてきた手法と同等、もしくはそれ以上の精度で対象問題に対する数値解を得られるレベルにあることが示された。
表(1).2 水平 2 層モデルを用いた問題設定における震源パラメータと物性値

6赤線：E-wave FEM
青線：差分法
点線：解析解図(1).1 水平二層モデルを用いた問題設定における原点での速度波形の比較

図(1).2 Moon Basin 問題における速度境界形状、観測点位置、及び震源パラメータ7

青線：真の境界形状

堆積層
母岩

E-wave FEM の計算用メッシュ

差分法の計算用格子点

図(1).3 y=0 平面における Moon Basin 問題の地盤構造に対して生成された E-wave FEM と差分法要の計算モデルの比較

図(1).4 図(1).2 で示される観測点上での計算速度波形の比較（黒点線：差分法。赤線： E-wave FEM）

8

観測点 B

観測点 C

観測点 D

時刻歴速度波形

振幅相対差=
疑似速度応答スペクトル|FDM-E-wave FEM|/E-wave FEM

図(1).5 観測点 B・C・D での速度波形、疑似速度応答スペクトル(h=0.05)とその振幅相対差（黒点線：差分法。赤線：E-wave FEM(右図のみ相対差)）9赤線： E-wave FEM（1 波長当たり 5 要素）
青点線： E-wave FEM（1 波長当たり 10 要素）

図(1).6 E-wave FEM においてメッシュサイズを変更して計算した観測点 D での速度波形の比較①-(c). 面的な広がりを持った震源モデルに対する地震動計算の精度検証
国の被害想定における長周期地震動評価では、面的な広がりをもった震源モデルに対する応答に対する計算が必要となる。面的な広がりをもった震源モデルを差分法、有限要素法などの数値計算に導入する方法としては、(i) 一定の間隔で分布する点震源の集合体を一つの震源モデルとみなす方法(ii) 節点分割法[5]の二つが考えられる。節点分割法は、断層面上に位置する節点の両側での仮想的な食い違い変位を考え、それと等価な体積力に対する応答を計算する方法である。有限要素法において、非構造要素により断層面位置
節点を配置し、節点分割法を導入することで、断層面上の連続的なすべりを点震源の集合を考える場合と比べてより精度よく考慮することができる。このようなアドバンテージは、震源位置が観測点と比較的近い場
合に発揮される。国の長周期地震動評価においても、節点分割法に基づいた震源モデルの導入を前提とした検討が進められている。E-wave FEM では、(i)(ii)の双方を導入可能となっているが、本節では、まず E-wave
FEM を用いて(i)のアプローチで従来の数値計算手法と同等の計算ができることを確認した上で、(ii)のアプ
ローチで得られる結果との差を比較する。これらの検討を通じて、C.「震源断層のモデリング」に関する精度検証を行う。
本節においては、面震源を入力する断層面形状の詳細な再現が重要なことから、比較対象の従来法として、大規模並列計算に実績があり、格子幅を細かくした計算に向いている差分法コードを用いることとした。その差分法コードの概要について述べる。ここでは[12]などに従い、食い違い格子を用いることとし、差分近似精度は、空間 4 次、時間 2 次とした。非弾性減衰の効果は、[13]の方法に基づいて考慮した。振幅の減衰性を考慮する参照周波数を f0、時間ステップを Δt、クオリティ・ファクターを Q(x,y,z)とすると、各時間ステップの速度と応力を計算する際に、次式で表される係数 A(x,y,z)を 1 つ前の時間ステップの値に乗じることで簡易に非弾性減衰の効果が評価される。

 − πf 0 ∆t A(x, y, z ) = exp  Q ( x, y , z ) (1-11)

10この方法では P 波の減衰 Qp と S 波の減衰 Qs を独立に与えることはできない。また、Q 値の周波数依存性
として𝑄𝑄 = 𝑄𝑄0 ⋅ 𝑓𝑓が仮定される。計算領域の側面と下面では，有限の領域を打ち切るために、[14]の one-way
波動方程式に基づく無反射境界条件と、[7]の吸収境界条件を併用した境界条件を設けた。[7]による吸収境

界は、吸収率 W を𝑊𝑊 = 𝑒𝑒𝑒𝑒𝑒𝑒(−𝛼𝛼 2 (𝐽𝐽0 − 𝑗𝑗)2)で表し、境界に近づくに従い指数的に W を増加させる。α は吸収率の大きさをコントロールする係数、J0 は吸収領域幅を表す。J0 を大きくするほど効果は大きく、対象とする計算問題に応じて J0 の幅を決める。[16]の方法に基づき、応力の各成分に対して、対応するモーメントテ
ンソルを入力することにより、食い違い震源断層を考慮する。

計算効率の向上のために、領域分割に基づいた並列計算を行えるようになっている。3 次元の計算領域を
NS、EW、UD の 3 方向に分割し、分割した小領域の演算を複数のノードに割り当てて計算負荷を分散させた。計算ステップ毎に MPI（Message Passing Interface）を用いた非同期通信を行い、各小領域の袖領域に設けたオーバーラップグリッドの速度と応力の値を、隣接するプロセス間で送受信して共有した（例えば[17][18]）。ノード内の並列計算もプロセス並列とした flat-MPI により行う。
面的な広がりを持った震源モデルに対する応答を計算するための問題設定は、[19]の N22 モデルに従うものとする。地盤モデルの概要を図(1).7 に、物性値を表(1).3 に示す。+X は北、+Y は東、+Z は鉛直下方を表す。地盤モデルの Q 値は∞であり、計算では Q=9999.9 と入力して、減衰の影響は考慮しない。断層長さ6km、断層幅 6km の低角逆断層の面震源を外力とした。断層の基準点は図(1).7 の①であり、走向 115deg.、
傾斜角 40deg.、すべり角 70deg.とした。破壊開始点を(0, 0, 6) (km)（モデル中央の深さ 6km の位置）とし、同心円状に破壊伝播速度 3km/s で破壊を伝播させた。すべり速度関数𝑆𝑆̇は次式を用いた。

 t − S ( X , Y=
, Z , t ) S0  2 ⋅ e T TtS ( X , Y , Z , t ) = 0（τ≧0）(1-12)(1-13)（τ＜0）

t は時刻(s)を表す。So と T は定数であり、それぞれ 1m、0.1s とした。
計算結果は、震央（0, 0, 0）と（0.6, 0.8, 0）を結ぶ自由地表面上の測線において、震央距離-10km～10kmの範囲を対象に、2km 間隔の計 11 点で出力した。[19]と同様に、X 成分と Y 成分として出力される計算結果を、観測点がすべて（0, 0, 0）と（0.6, 0.8, 0）を結ぶ直線長にあることを踏まえ、次式で Radial 成分（R）と Transverse 成分（T）に変換した。

R=−0.6 X − 0.8Y
T=
+0.8 X − 0.6Y

（X＜0，Y＜0 の地点）

(1-14)

R=+0.6 X + 0.8Y
T=
−0.8 X + 0.6Y

（X≧0，Y≧0 の地点）

(1-15)

上下成分は Up を＋、Down を－となるように符号を変えた。計算結果に対して図(1).8 示す周波数特性のローパスフィルタ処理を施した。
差分法計算において、計算領域は、境界処理に必要な部分を除いて（X, Y, Z）=（30, 30, 17）(km)とした。
計算領域の側面と下面に 60grid 幅（水平 3km 幅, 鉛直 1.2km 幅）の吸収領域を設けた。先述したように、
面震源を入力する断層面形状の詳細な再現を目的とし、X 方向と Y 方向の格子間隔は 50m、Z 方向では 20m
11とした。時間刻みは 0.001 秒とした。面震源は、点震源の重ね合わせにより考慮し、50m メッシュで分割した 14,641 個の点震源に対してダブルカップル力を与えた。X, Y, Z の 3 方向をそれぞれ 4 分割した 64 並列により行った（例えば[17][18]）。SX-Aurora TSUBASA（NEC 社製）を計算に用いた。E-wave FEM にお
いては、解析領域は、水平方向を 30km×30km、6 Hz の波に対して、1 波長当たり最低 5 要素が確保される
ようにメッシュ生成した。計算時間は 0.02 秒刻みで 10 秒間とした。総ステップ数は 500 である。計算はOakforest-PACS を用い、8 ノード 64MPI プロセスの並列数で実行した。E-wave FEM では、震源モデル
を(i)「一定の間隔で分布する点震源の集合体を一つの震源モデルとみなす方法」と(ii)「節点分割法」の二種類の方法で導入する。前者を Case 1、後者を Case 2 として両方の結果を示す。
E-wave FEM（Case1）と差分法の比較として、速度波形を図(1).9、フーリエスペクトルを図(1).10 に示す。E-wave FEM と差分法の結果は全ての成分・地点でよく一致している。[19]においても点震源を用いたアプローチがとられており、これらの計算結果の数値データとは比較できなかったものの、波形を重ねるとよく一致した。ELAPSE TIME は 02:14:11 であった。次に、E-wave FEM の Case1 と Case2 の速度波形の比較（図(1).11）においては、E-wave FEM と差分法の位相は Radial 成分と Vertical 成分では概ね一致しているが、Transverse 成分では違いが見られる。最大値については大半の地点・成分で Case2 の方が小さい。Case1 と Case2 のフーリエスペクトルの比較を図(1).12 に示す。Radial 成分と Vertical 成分については、Case2 は Case1 に概ね一致する。詳細に見ると多少の違いがあり、例えば、Radial 成分の+000 地点
から-004 地点のピーク周波数付近では Case2 が Case1 よりも小さい。Transverse 成分については、+010
地点では Case2 は Case1 に比較的類似しているが、その他の地点では Case1 と Case2 の違いが Radial 成
分と Vertical 成分よりも顕著である。ELAPSE TIME は 02:19:44 であった。
まとめると、E-wave FEM で(i)点震源の集まりを用いて震源モデルを表現した場合は、全ての成分で差分法によく一致する結果が得られた。E-wave FEM で(ii)節点分割法を用いた場合、Radial 成分と Vertical 成分については差分法と概ね一致する結果が得られたが、Transverse 成分については差分法との違いが確認された。点震源の集合と節点分割法を用いた場合の結果の差は無視できない場合があり、面的な広がりをもった震源をより忠実に導入できる節点分割法の有用性が示唆される一方で、結果の差の原因について今後精査が必要と考えられる。

図(1).7 地盤モデルの概要（1/4 対象部分）（[19]）

12

表(1).3 地盤モデルの物性値（[19]）図(1).8 ローパスフィルタの周波数特性（[19]）13図(1).9 E-wave FEM（Case1：黄破線）と差分法（青実線）の速度波形（Radial 成分、Transverse 成分、Vertical 成分）

14図(1).10-1 E-wave FEM（Case1：黄破線）と差分法（青実線）のフーリエスペクトルの比較（-010 地点～+004 地点）

15図(1).10-2 E-wave FEM（Case1：黄破線）と差分法（青実線）のフーリエスペクトルの比較（+006 地点～+010 地点）

16図(1).11 E-wave FEM の Case1（黄破線）と Case2（青実線）の速度波形（Radial 成分、Transverse 成分、Vertical 成分）

17図(1).12-1 E-wave FEM の Case1（黄破線）と Case2（青実線）のフーリエスペクトルの比較（-010 地点～+004 地点）

18図(1).12-2 E-wave FEM の Case1（黄破線）と Case2（青実線）のフーリエスペクトルの比較（+006 地点～+010 地点）①-(d). 表面波が長時間持続する水平 2 層問題における解析解との比較
我々の目的は長周期地震動の評価なので、長く持続する表面波が主体となるような問題を長時間安定して解ける必要がある。このような問題においては、計算領域端部における反射波の処理精度、媒質における減衰の計算精度が波形計算の精度に大きな影響を及ぼしうるため、D.「モデル端部における反射波の処理」とE.「減衰のモデリング」に関する精度検証が必要となる。ここでは、D・E についての精度検証と、その結果必要となった D についての機能追加について述べる。
まず、解析解の存在する水平 2 層問題で、地表面波形において長く持続する表面波が主体となるような問題を設定し、解析解と E-wave FEM による有限要素解を比較する。堆積層が地表面から厚さ 1.5 ㎞、それ以深は母岩という設定とする。物性値と震源パラメータを表(1).4 に示す。先に扱った水平 2 層問題と比較し、2 層の物性値のコントラストが大きい、1 層目の中に震源を設置している、などの違いがある。このような設定下では、堆積層と母岩の速度境界面において振幅の大きい波が繰り返し反射し、複雑かつ長時間持続 す る 波 動場 が 形 成され る 。E-wave FEM での 計 算 に おい て は 、計算 領 域 を -80km<=x<=48km 、 32km<=y<=32km、-56km<=z<=0km と設定し、これまでと同様に 1 波長当たり最低 5 要素という計算条
件で計算する。この問題設定により得られる地表面波形の例として、原点での速度波形を図(1).13 に示す。
解析解を示す黒点線によって示されるように、図(1).1 の波形に比べ、波がはるかに長い時間持続していることがわかる。E-wave FEM の計算波形との比較においては、後続波において解析解との波形のずれが目立っていることがわかる。その後の何度かの試行計算により、波形のずれ方が有限要素モデルのモデル領域の変化に伴い変化することが分かったため、領域端部における反射波が E-wave FEM にもともと搭載される粘性境界条件[6]（表(1).1D）だけでは処理しきれていないことが示唆された。そこで、領域端部に吸収帯、あ19

るいはスポンジ領域と呼ばれるバッファ領域を設け、この領域で弾性波の振幅を人工的に抑制する手法[7]を導入することとした。この手法では、(𝑖𝑖)(1-16)(𝑖𝑖)

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

(𝑖𝑖)

(𝑖𝑖)

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

𝑊𝑊 = �

exp�−𝛼𝛼(𝐽𝐽0 − 𝑗𝑗)� (𝑗𝑗 ≤ 𝐽𝐽0 )1 (𝑜𝑜𝑜𝑜ℎ𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒)(1-17)(𝑖𝑖)

𝑗𝑗 = 𝑑𝑑 /𝑑𝑑0のように定義される。ここで、𝐽𝐽0 , 𝛼𝛼, 𝑑𝑑 (𝑖𝑖) , 𝑑𝑑0はそれぞれ、吸収帯の厚さに関わるパラメータ、吸収帯内のダン

ピングの空間変化率を決めるパラメータ、i 番目の節点の最近傍モデル境界からの距離、正規化用距離である。𝑑𝑑0 としては水平方向の最大メッシュサイズを選択する。このようにして、領域端部での反射波が主たる計算領域に影響を与えないようにする。これは速度を未知数とする差分法スキームにおいて広く使用されているが、変位を未知数として解く有限要素法においても使われる[21]。E-wave FEM に𝐽𝐽0 ＝100, 𝛼𝛼 = 0.002というパラメータに基づき吸収帯を導入して計算した波形と、解析解による原点での波形の比較を、図(1).14aに示す。吸収帯を導入したことで、後続波において波形が解析解とよく一致するようになった。x 成分の計
算速度波形について、𝐶𝐶 ∗スコアは 8.81 から 9.34 に改善した。80 秒前後でやや目立つずれは、1 波長当たり

最低 10 要素として計算したところほとんど目立たなくなり（図(1).14b）、𝐶𝐶 ∗スコアも 9.68 とさらに改善したため、離散化の影響によるものと考えられる。これより、吸収帯の導入により D.「モデル端部における反射波の処理」が上手くいくようになることが確認できたといえる。

続いて、同様の問題でレイリー減衰による E.「減衰のモデリング」の精度を確かめる。ここまで示した例
においては、Q 値が 1000 と大きかったため、減衰が計算波形に及ぼす影響は小さかった。そこで、新たに堆積層と母岩の Q 値を P・S 波共通でそれぞれ 160、400 として、他の設定は変えずに波形計算を行った。
図(1).15 に、1 波長当たり最低 10 要素とした場合の原点での波形の比較を示す。Q 値が小さい場合でも Ewave FEM の計算波形と解析解がよく一致（x 成分の計算速度波形 は𝐶𝐶 ∗スコア 9.70）していることから、レイリー減衰による E.「減衰のモデリング」が、表面波の伝播が長時間にわたる場合の計算においてもうまく機能していることが示された。

20

表(1).4 長く持続する表面波が発生するような水平二層問題における震源パラメータと物性値

図(1).13 粘性境界条件のみを用いた場合の原点における速度波形（黒点線：解析解、赤線：E-wave FEM）

21

a

b

図(1).14 吸収帯を用いた場合の原点における速度波形
a: 1 波長当たり最低 5 要素とした場合、b:10 要素とした場合（黒点線：解析解、赤線：E-wave FEM）

22

図(1).15 吸収帯を用い Q 値を小さくした場合の原点における速度波形（黒点線：解析解、赤線：E-wave FEM）①-(e). 2013 年 4 月 13 日に淡路島で発生した地震に対する地震動計算
ここまで、理想的・単純な問題設定において、A-E の 5 項目について確認、E-wave FEM の精度検証・機能追加を進めてきた。続いて、より現実的な地盤モデルを対象にした長周期地震動計算への適用を行う。ここで検討すべきなのは、・E-wave FEM を現実的な地盤モデルにおける計算に適用した場合の課題の抽出と解決である。結果として B.「計算モデルの生成」について課題が抽出されたため、機能追加の対応したことについて説明する。
ここで対象とする問題は、近畿地方で過去に発生した地震の再現計算である。実際に推定されている震源パラメータに対する地震動を計算するため「再現計算」としているが、ここでは現実的な地盤モデルにおける計算に適用した場合の課題の抽出を目的とする。対象とする地震は 2013 年 4 月 13 日に発生した淡路島付近の地震(Mj 6.3）である。防災科学技術研究所の F-net によるメカニズム解に準じて設定する点震源のパラメータと、点震源の位置と今回注目する観測点の位置を図(1).16 に示す。この周辺領域に対し、同研究所の J-SHIS データベースをもとに、震源を含む東西方向 180km、南北方向 150km、鉛直方向 50km の領域に地震波速度が一様な 33 の層からなる三次元不均質速度構造を設定することとした。対象領域において、
J-SHIS の地盤データにおいては、中央構造線付近に地震波速度境界の勾配がかなり大きい部分が存在している。E-wave FEM に搭載のメッシュ生成手法は背景構造格子を活用した四面体二次要素のメッシュ生成
アルゴリズムを採用しているが、本手法を標準的に用いたところ、共役勾配法を用いた連立一次方程式の求解において、反復の収束に失敗した。これはメッシュに非整合が存在する場合に見られる症状の一つであり、勾配が大きい速度境界を含む地盤データに対して整合的なメッシュを生成することに失敗したことがわかった。
勾配の大きい速度境界を含む地盤データに対して整合的なメッシュを生成するために、拡張機能の導入と機能追加を行った。図(1).17 にその概要を示す。図(1).17a は、地表面・速度境界に急勾配を持つ 2 層からなる地盤構造 DEM データの 2 次元模式図である。E-wave FEM の標準的なメッシュ生成手法では、DEM データを含むメッシュ生成対象領域に、立方体（図では正方形）からなる背景構造格子を生成し、各々の立方体内に地表面・速度境界が含まれるかを判定し、含まれる場合には地表面・速度境界に沿うように立方体を
23

四面体要素により分割する。しかし、図(1).17b の赤枠部分のように、地表面・速度境界が縦方向に複数の構造格子をまたぐ場合、この判定が正しくできなくなり、メッシュが上手く生成できなくなる場合がある。対処法として、図(1).17c のように、DEM データの底面からの高さをいったん圧縮してメッシュを生成し（この圧縮率を Squash factor と呼ぶ）
、再び元の大きさに引き延ばすことで、図(1).17d のように、元のメッシ
ュ生成手法を活用したまま、急勾配を含む DEM データに対してもメッシュ生成が可能となる。新たに発生する問題として、1 波長当たりの要素数のコントロールが難しくなることがあげられる。その理由は、標準的なメッシュ生成手法においては、背景構造格子の立方体サイズを制御することで 1 波長当たりの要素数をコントロールしていたが、Squash factor の導入により水平と鉛直方向で辺の長さが異なる構造格子がベースとなるからである。そこで、メッシュ生成後、希望する 1 波長当たりの要素数を満たさない部分の要素については、1 辺の長さが 1/2 の 8 個の四面体要素に再分割することとした。その際、周囲の要素も図(1).17eのパターンにしたがって再分割することで、メッシュ全体の整合性が保たれる。近畿地方の地盤構造に対して、Squash factor を 3 として要素再分割を伴うメッシュ生成により、0.25Hz までの地震波を対象に、1 波長当たり 5 要素を保証してメッシュを生成した。このようにして生成された計算モデルの例の概観を図(1).18 に示す。計算を行ったところ、期待される反復回数で連立一次方程式の求解が可能となったことを確認できた。これより、勾配が大きい速度境界を含む DEM データに対しても、1 波長当たりの必要要素数を保証して堅牢にメッシュを生成することが可能となったことを確認した。
現実的な地下構造モデルでの地震動計算の精度検証のため、E-wave FEM での数値解の適切な設定について検討する。まず、1 波長当たり 5 要素という設定が十分かの確認を行う。1 波長当たり 5 要素を保証してメッシュにおける計算においては、点震源の深さが 16.4 ㎞となり、設定したい値から 0.6 ㎞ずれることが分かった。これは、squash factor＝3 としたことで、縦方向の要素サイズが大きくなったことが原因と考えられる。1 波長当たり 10 要素を保証してメッシュを生成したところ、点震源を 17 ㎞の深さに設定することが可能であることが分かった。これにより、1 波長当たりの要素数を変更してメッシュサイズを変えると、点震源位置にもずれが出てしまうことが分かったので、
1. 1 波長当たり 5 要素＋震源深さ 16.4 ㎞と 1 波長当たり 10 要素＋震源深さ 16.4 ㎞
2. 1 波長当たり 10 要素＋震源深さ 16.4 ㎞と 1 波長当たり 10 要素＋震源深さ 17 ㎞の二段階に分けた比較により、メッシュサイズと震源深さの影響を切り分けて議論する。まず、1 の点震源深さをそろえてメッシュサイズを変更した場合の結果を図(1).19 に示す。
それぞれ 1 波長当たり 5 要素と 10要素の計算結果を表す黄線と緑線の比較において、特に右列に示される 0.25Hz のローパスフィルタをかけた場合だと、多少の振幅の違いがあるものの、黄線が重なって目立たないほどによく一致している。つまり
1 波長当たりの要素数の違いに関わらず計算速度波形がよく一致しており、例えば HYG018 点 xyz 三成分の速度波形における𝐶𝐶 ∗スコアはそれぞれ 9.99、10.00、10.00 となっている。2 での点震源深さに違いがある

場合の結果を図(1).20 に示す。この場合は、それぞれ点震源深さが 17km と 16.4km の計算結果を表す黄線と緑線の比較において、左列のフィルタ―をかけない場合と右列のフィルターをかける場合の双方で、振幅はよく一致しているものの、振幅のピークの時刻のずれが比較 1 と比べて目立っている。つまり震源深さの違いが位相の違いとなってかなり明確に表れているといえる。HYG018 点 xyz 三成分の速度波形における𝐶𝐶 ∗スコアはそれぞれ 9.87、9.99、9.95 と、比較 1 の場合と比べて小さい。これらの結果より、1 波長当たり5 要素という設定は、メッシュサイズという観点からは既往手法との比較に十分と思われるものの、点震源位置の変化による影響が無視できないため、今回は 1 波長当たり 10 要素という設定により得られた数値解を最終的な結果と考える。1 波長当たりの要素数を変えずに点震源位置を想定する位置に近づけることは、
今回導入した要素再分割をその位置の周辺部分に適用することで可能となる。ここまでの計算では、モデル端部における反射波の処理のために、E-wave FEM に従来組み込まれている粘性境界のみを用いたが、GMSとの比較の際には吸収帯も導入することでより反射波処理の性能を向上させることを試みる。図(1).21 に、
24

吸収帯のパラメータを𝐽𝐽0 ＝120, 𝛼𝛼 = 0.001として計算した場合と、吸収帯を導入しない場合の計算結果の比較を示す。両ケースともメッシュサイズは 1 波長当たり 10 要素と変わらないが、吸収帯を導入した場合では吸収帯の中に震源が含まれないよう、鉛直方向のモデル厚さを大きくしている。それぞれ吸収帯を使用する
場合としない場合の計算結果を表す黄線と緑線の比較において、左列のフィルタ―をかけない場合と右列のフィルターをかける場合の双方で、反射波の影響がない直達波の部分の波形はおおむね一致しているが、反射波の影響が表れうる後続波において、吸収帯を導入したことによる波形の変化がみられる。吸収帯を導入した計算結果の方が、反射波の影響を低減したより精度の良い数値解となることが期待されるが、それを直接的に検証することはまだできておらず、今後の課題となっている。1 波長当たり 5 要素のケース、10 要素のケース、10 要素で吸収帯を導入したケースでは、モデル自由度はそれぞれ約 4 億 3 千万、33 億 2 千万、64 億 7 千万となった。すべてのケースにおいて時間ステップ幅 0.01 秒で 6,000 ステップ分計算した。用いた計算資源はそれぞれ Oakforest-PACS の 128 ノード×19 分、128 ノード×3 時間 30 分、128 ノード×8 時間となった。
A-E の 5 項目に関するここまでの精度検証と機能追加により、E-wave FEM の国の長周期地震動評価のための計算に使えるツール（プリプロセス＆計算コード）としての準備ができたことを確認した。

図(1).16 2013 年 4 月 13 日に発生した淡路島付近の地震(Mj 6.3）の震源パラメータ、及び地図上の震源と観測点位置
黄色点で表される対象領域内に K-NET・KiK-NET 観測点のうち、黒点のものに注目する。

25

a

b

c

d

e

図(1).17 メッシュ生成手法への機能追加の概要
a. 地表面・速度境界に急勾配を持つ 2 層からなる地盤構造 DEM データの 2 次元模式図。b. Ewave FEM の標準的なメッシュ生成手法で用いられる背景構造格子。赤枠の部分などで、速度境界が縦方向に複数の構造格子をまたぐ形になっている。c. DEM データの底面からの高さを1/2 としてから、つまり Squash factor=2 としてメッシュ生成する様子。d. c で作られたメッシ
ュを、元の DEM と整合するように再び引き延ばしたもの。e. 二次元要素の再分割パターン。

大阪平野

淡路島

図(1).18 近畿地方に対して生成された有限要素メッシュの概観
青に近いほど浅い層、赤に近いほど深い層を表している。

26

図(1).19 1 波長当たり 5 要素＋震源深さ 16.4 ㎞と 1 波長当たり 10 要素＋震源深さ 16.4 ㎞の計算波形比較27

図(1).20 1 波長当たり 10 要素＋震源深さ 16.4 ㎞と 1 波長当たり 10 要素＋震源深さ 17 ㎞の計算波形比較28

図(1).21 吸収帯のパラメータを𝐽𝐽0 ＝120, 𝛼𝛼 = 0.001として計算した場合と、吸収帯を導入しない場合の計算波形比較。

①-(f). 相模トラフ地震に対する関東地方での長周期地震動計算への適用
ここまでの検討で、E-wave FEM が検討項目 A-E に関する精度検証・機能追加を通じ、理想的・現実的な地盤構造モデルにおける長周期地震動計算を計算するのに十分な機能を有することを確認した。最後に、
E-wave FEM を「富岳」上にて用い、相模トラフ域を震源とする関東地方における長周期地震動のシミュレーションに適用する。同じタイプの地震として元禄地震や大正関東地震が知られており、国の被害想定において重要な検討対象となっている。実際に国の評価で用いる地盤構造データと同じものでメッシュ構築を行い、必要な計算設定を確認の上、「富岳」上で長周期地震動計算を実施する。本計算を通じて、E-wave FEM
29を用いて国の被害想定の算出に対して必要な規模の地震動計算を実施できるようになっていることを確認する。
対象とするのは、伊豆半島から北関東、房総半島までを含む領域である。この領域に対し、防災科学技術研究所の J-SHIS データベースをもとに、地震波速度が一様な 33 の層からなる三次元不均質速度構造を設定する。大陸マントルとフィリピン海プレートの海洋性地殻の速度境界をプレート境界と考え、そこに震源断層を設定する。採用した震源断層モデルでは、神奈川県南部、東京湾入口付近、房総半島南端の地下に強震動生成領域が設定されており、空間的な広がりを持ち、破壊伝播を伴う。震源断層モデルは、C.「震源断層のモデリング」にて従来法の結果の差についての確認を行った節点分割法[5]により、有限要素方程式の外力ベクトルへと変換して入力する。震源時間関数は、任意の震源時間関数に対する応答を合成できるように、三角形型のものを用いることとした。震源断層モデルは検討中のものであり、時間方向にはインパルス応答的な計算となっているため、計算された地震動は直接的には地震被害に対する情報を持たない。0.5Hz までの弾性波を対象に、1 波長当たり 5 要素の計算条件で、有限要素メッシュを生成することとした。この計算
条件は、①-(e)の現実的な地盤構造における検討で、1 波長当たり 5 要素と 10 要素の波形の一致度が高かった（図(1).19）ことにより採用することとした。時間ステップ幅を 0.02 秒とし、104 ステップ、つまり 200
秒分の地震動を計算することとした。
吸収帯の設定については、モデル領域や震源位置との関係で、問題によって適切な設定が変わるため、今回の問題設定について必要な設定を検討した。検討開始時のモデルは、平面直角座標 9 系において、東方向を x 方向、北方向を y 方向として、原点位置を(x,y)=(-170km,-160km)の位置にとり、モデル領域を東西方
向 294 ㎞、南北方向 336 ㎞、深さ 50 ㎞としたものであった。モデル領域に変化を与えた場合に、波形にどの程度変化がみられるかを調べることにより、反射波が計算波形にどの程度影響するかを調べることができる。検討開始時のモデルと、このモデル領域を東西南北の水平方向各モデル端部を法線方向に 4.5 ㎞分小さ
くした場合の速度波形の比較の例として、房総半島東部のある観測点での結果を図(1).22a に示す。反射波が原因と思われる波形のずれが多少目立っており、三成分の𝐶𝐶 ∗スコアはそれぞれ 9.92、9.96、9.83 となった。同じメッシュ生成条件の下で計算していることを考慮すると若干低めの値となった。その他の点においても、場所によって若干ずれが目立つ観測点が存在した。吸収帯の有無・モデル領域の双方について検討を重ねた結果、吸収帯については𝐽𝐽0 ＝120, 𝛼𝛼 = 0.001（本モデルにおいては水平方向厚さ 33.6 ㎞、鉛直方向厚

さ 67.2 ㎞）のパラメータにより設置するとともに、設定した震源モデルに対する応答が主たる計算領域に影響を与えないように、
モデルの南端を y の負の方向に 33.6 ㎞伸長したモデル領域を設定することとした。

このモデルにおける計算速度波形と、モデル領域に同様に 4.5 ㎞の変化を与えた場合の速度波形の比較の例
を図(1).22b に示す。見た目にも波形のずれが目立たなくなっており、三成分の𝐶𝐶 ∗スコアはそれぞれ 9.99、
9.99、9.99 と改善した。最終的なモデルの自由度は 588 億程度となった。計算モデルの概観を図(1).23 に示す。このモデルに対し、「富岳」4,096 ノードを用いて約 5 時間で計算を終了した。最初の 10 ステップにおいて、共役勾配ソルバー部分における対ピーク比計算性能は 7.8％程度であった。図(1).24 に計算された波動場（変位）を可視化したものを示す。先に示した波形と合わせ、設定した震源モデルに対し、不自然な成分や反射のない波動場が計算できていることが確認できる。以上の検討から、応答を計算すべき震源断層モ
デルが策定されれば、「富岳」を用いた関東地方における長周期地震動の応答を計算することができる状況となったことを確認した。次年度、国の担当機関において「富岳」の計算資源を確保して長周期地震動評価
が実際に進められる予定となっているが、このような計算例を踏まえ、その際に必要な計算資源の見積への協力も進めてきている。このように、国の被害想定の検討に「富岳」にて E-wave FEM が使用される準備がほぼ整ったといえる。同規模の計算を第二階層の計算機で行うことは現実的ではない（例えば OakforestPACS の場合、システム全体の約半分を 10 日程度占有しなければならないことになる）。またこの状況は、
E-wave FEM を用い、今後 2.5Hz 程度までといったさらに高周波までの地震動計算を行うことで、非線形30地盤増幅計算プログラム「STRIKE」を用いたローカルな浅層地盤における非線形地盤増幅計算のための入力地震動を作成する準備が進んでいることにも相当する。
今後必要な検討事項としては、節点分割法による震源断層の入力方法が、点震源の集合を用いる方法とどのような関係にあるかということがある。これまでの国の被害想定における長周期地震動評価においては、
後者の方法が採用されているため、これらの方法の違いが計算結果にどのような差を及ぼすのか精査を行う
予定である。

図(1).22 a. 検討開始時のモデル（赤線）と、このモデル領域を東西南北の水平方向各モデル
端部を法線方向にそれぞれ 4.5 ㎞分変化させたモデル（黒線）での、房総半島東部のある観測点での速度波形の比較。b.最終的に検討に用いることにしたモデル（赤線）と、このモデル領域を東西南北の水平方向各モデル端部を法線方向にそれぞれ 4.5 ㎞分変化させたモデル（黒線）での、同様の比較。
31

a

b

図(1).23 関東地方の 33 層からなる不均質地盤構造に対して生成した計算モデル
a. モデルの概観。側面において、色が赤に近いほど地震波速度が小さく、青に近いほど大きい。
b. a 中赤線部での断面図において可視化された有限要素メッシュ。

富士山

房総半島

図(1).24 計算された波動場（変位）の可視化図
変位が地表面の変形として表されており、さらに振幅が大きい部分は赤色で強調されている。

➁ E-cycle FEM を用いた弾性グリーン関数ライブラリの計算と公開準備
断層すべりに対する弾性グリーン関数は、地震時の断層すべり推定や地震間の断層固着度推定に加え、津波シナリオ計算など防災に関しても様々の用途があり、沈み込み帯地震に関する地震学・地震工学研究における重要なコンポーネントと位置付けられる。線形性により、一度計算しておけば重ね合わせによりいろい
ろな現象の解析に適用できる。[21]などによる計算コードの整備により、半無限弾性体での弾性グリーン関
数計算が広く行われ、測地観測と組み合わせた断層すべり逆解析などの分野で成果を挙げてきた。一方で、
地表面形状はもとより、地下の弾性構造についても、地震探査やトモグラフィにより情報が蓄積されるようになってきた。また地上の GNSS 地殻変動観測が常時行われるようになったことに加え、海底地殻変動観測や干渉合成開口レーダー (InSAR) による地殻変動解析がより多くの地域で可能となってきた。これらの観測から得られるデータがカバーする広い範囲での地殻変動データを整合的に扱おうとする場合、地形、弾性パラメータ、地球の楕円度を考慮したグリーン関数の計算、および断層すべり逆解析、津波計算など望まし
いと考えられるようになってきた。
32

グリーン関数の使用に関しては、一度計算をしておけば後はユーザが簡便な計算で活用できるという特徴を考えれば、計算は各研究者が各々で行う必要はなく、多くの人に使われるために十分な要求を備えたグリーン関数のライブラリが公開されていればよい、と考えることもできる。そのため、課題参加者らは重点課
題時代から、地殻変動計算プログラム「E-cycle FEM」を用い、日本の主要な沈み込み帯に対し、3 次元弾性不均質性を考慮した弾性グリーン関数の有限要素計算を行い、ライブラリとして公開するための検討を進めてきた。本課題においては、重点課題時代のプログラム開発・手法検討・計算のフェーズから、計算・データ公開・活用のフェーズに入っている。本項目では、E-cycle FEM とグリーン関数計算の仕様の概要について述べた後、グリーン関数ライブラリのデータ公開・活用の現状について紹介する。
➁-(a). E-cycle FEM の概要
E-cycle FEM は、断層における食い違いに対する静的な弾性応答・準静的な粘弾性応答を計算する有限要
素計算プログラムである。ここでは簡単のために今回のグリーン関数ライブラリの作成に使用する静弾性応
答のみについて説明する。静弾性の支配方程式に対し有限要素定式化を行い、境界条件を考慮して整理すると、問題は以下の一次連立方程式に帰着する：(1-18)𝐊𝐊𝐊𝐊 = 𝐟𝐟

ここで K、u、f はそれぞれ全体剛性マトリクス、変位ベクトル、外力ベクトルである。断層におけるすべり
量を節点分割法[5]により f に変換する。E-wave FEM と同様、本一次方程式の求解が最も計算資源を必要とする部分となる。求解を高速に行うため、高性能なソルバーが導入されている[22]。ベースとなるアルゴ
リズムが E-wave FEM と共通であり、同様に良好なスケーラビリティが得られることに加え、本式では質量項が存在しない分反復法の収束が比較的悪いことを考慮し、これに加え代数的マルチグリッド法も併せて用いることで、収束性のさらなる改善をはかっている。このような工夫により、1012 オーダーの自由度の計算モデルに対しても、現実的な時間での求解が可能となる[22]。固体地球分野においては 107 オーダーの自由度の静弾性応答計算（例えば[23]）が一般的であり、本手法において可能な計算の詳細さは他に類を見ない。
➁-(b). グリーン関数ライブラリの概要
ここでいうグリーン関数とは、対象とするプレート境界を離散化する小断層に対する弾性応答である。グリーン関数の計算は南海トラフと日本海溝を対象として計算を行った。日本列島においては、統一的な地震
波速度構造モデルとして「全国一次地下構造モデル（暫定版）」[24][25]が構築されている。地理座標系で定義されるこれらのデータをデカルト座標に変換して、E-wave FEM と共通のメッシュ生成手法を用いて南海トラフ域（図(1).25a）と日本海溝域（図(1).25b）のメッシュを生成した。メッシュ生成のために用いたパラメータは、日本海溝域での M9 地震に対する粘弾性応答を 1%以下の誤差レベル（応答レベルが小さな部分については期待される観測精度以下の誤差レベル）で計算するための設定を調べた[26]に基づいて設定した。これにより、層境界の水平方向要素サイズを 1 km、震源域を中心として領域を 2500 km×2500 km×1100 km としてメッシュを生成した。ただし鉛直方向については、対応する要素サイズが 1.5 ㎞となっている部分がある。これは、使用している背景構造格子ベースのメッシュ生成手法[2]において、一部急勾配を含む DEM データに対しメッシュを整合的に生成するため、鉛直方向が水平方向に比べて 1.5 倍の長さであるボクセルをベースとすることが必要となったからである。また、最浅層の S 波速度が 3550m/s であった[26]のに比べ、今回使用した地震波速度構造においてはそれが大幅に小さい 700－2000m/s 程度の層が存在した
ため、水平方向要素サイズが 1 ㎞で十分であるかをメッシュサイズ ds を変えた事前解析により検証した。
33

具体的には、上盤側において S 波速度が 700－2000m/s 程度と小さい堆積層より浅い層における水平方向要
素サイズを再分割により 0.5km、0.25 ㎞と変化させたメッシュを準備し、[26]で用いたものと同様の断層すべりを日本海溝域モデルに入力し、海域の観測点における弾性応答の変化が 1％の範囲に収まる要素サイズを採用することとした。この結果、上盤側の堆積層より浅い層の水平方向メッシュサイズのみ 0.5 km（鉛直方向については同様に一部 0.75 ㎞）とした。計算モデルの自由度と要素数は、南海トラフ域モデルはそれ
ぞれ約 59 億、約 15 億、日本海溝域モデルはそれぞれ約 56 億、約 14 億となった。小断層として、3 次のBi-cubic B スプライン形状を持った単位すべりをプレート境界面上に配置した構造格子上の各点に入力する。2 方向のすべり方向を考慮する。1 方向は格子点におけるプレート境界面への接面に沈み込み方向を投影した方向とし、もう 1 方向は格子点におけるプレート境界面に対する法線ベクトルと沈み込みベクトルの双方に直行する方向とする。西南日本の南海トラフ域におけるプレート境界の主だった領域に Bi-cubic B ス
プライン形状の小断層を導入するための格子点を配置すると、その個数は 460 個となり、沈み込み方向とその垂直方向とを合わせて 920 個のグリーン関数を計算することとなった（図(1).25c）。東北日本の日本海溝域に対しても同様に格子点を配置すると、その個数は 405 個となり、沈み込み方向とその垂直方向とを合わ
せて 810 個のグリーン関数を計算することとなった（図(1).25d）。これらの枠組みは重点課題時代に構築済みであったが、もともと堆積層の水平方向メッシュサイズを 1 ㎞としていたのを、1.5 ㎞に再分割すると判断した後の再計算を、南海トラフ域についてだけ本課題において実施した。計算には Oakforest-PACS の
512 ノードを用いて 16 時間余りかかった。計算内容については、査読付き国際誌に報告済みである[27]。
➁-(c). グリーン関数ライブラリの公開と活用
グリーン関数ライブラリのデータ公開に際しては、計算したグリーン関数本体（1 km×1 km グリッドの
地表面変位・傾斜・体積ひずみ）に加え、● グリッド上のデータから任意の観測点でのデータに補間するプログラム● グリーン関数を重ね合わせて任意の断層すべり分布に対する地表面応答を計算するプログラムを同梱した環境一式を用意している。公開用の WEB サイトは https://www.jamstec.go.jp/feat/gflsz/に開設済みであるが、2021 年 5 月現在、所属機関の事情によりデータセット本体の公開自体はできていない状況にある。状況が改善し次第速やかに公開する予定である。
データ公開に先立って、共同研究という形で本課題外の研究者によるグリーン関数ライブラリの活用が進められている。一例として、南海トラフ域における地震間（非地震時）の陸・海域の測地観測データを用いたすべり遅れレートの推定への適用例[28]を図(1).26 に示す。ここでは、図中 a・c 黒矢印および b・d 黒太線で表される地表面変位レートの測地観測データから、色コンターで表されるプレート境界でのすべり遅れ
レート分布の推定が行われている。変位レートデータとすべり遅れレートパラメータの線形関係を表現する
行列の作成にグリーン関数ライブラリが使用されている。プレート境界断層に対し、3 次元不均質地下構造の影響を考慮した地殻変動を手軽に計算できる本ツールは、特にこのような陸・海域双方の応答変位を考慮する必要がある場合に需要があると考えられ、本例以外にもすでに複数の機関と連携した研究を進めている
ところである。
➁-(d). 終わりに本グリーン関数ライブラリは、固体地球科学分野での知見に基づき、本分野での地殻変動・津波波源に関する研究に耐えうる計算精度を担保することを想定して計算されている。本研究で構築した枠組みをもとに、
今後は津波防災に資するような信頼性を担保したグリーン関数計算を行うことも重要となってくる。そのためには、グリーン関数計算を行う際の地下構造モデル、メッシュサイズ、領域サイズ（モデル端部の境界条
件の影響）などを再検討する必要があると考えられる。本課題においては、静弾性・粘弾性応答関数の参照
34

解創出を目標の一つに掲げており、令和 3 年度以降に取り組む予定となっている。その理由は、参照解創出のためには、より詳細な地下構造に対し超高精細なメッシュと超大規模領域を採用した計算を行う必要があり、「富岳」でしかできない計算内容となると考えられるからである。
a

b

c

d

図(1).25 有限要素モデルの概観
a は南海トラフ域、b は日本海溝域を表す。赤線は地下構造（P 波速度）のプロファイルの可視化を実施した場所を表すが、その内容については本稿では割愛する。cd. グリーン関数を計算する対象領域の概観。c は南海トラフ域、d は日本海溝域を表す。灰色の点は各小断層の中心点を表す。三角形と逆三角形はそれぞれ GEONET の観測点と海底地殻変動観測点のうち代表的なものを示す。Pc は有限要素モデルの中心点を表す。

35

a

c

d

b

図(1).26 グリーン関数ライブラリの南海トラフ域における地震間（非地震時）のすべり遅れ
レート推定への適用例[28]。a. よく使われる推定手法（すべり遅れ分布の滑らかさに対する事前情報を加えたもの）による推定結果（色コンター）、観測水平変位レート、計算水平変位
レート（矢印）。b. その推定誤差（色コンター）、観測鉛直変位レート、計算鉛直変位レート（縦棒）。c. 異なる推定手法（すべり遅れ分布のパラメータ数自体も未知数とするもの）による推定結果（色コンター）、観測変位レートと計算変位レート（矢印）。d. その推定誤差（色コンター）、観測鉛直変位レート、計算鉛直変位レート（縦棒）。

参考文献[1] Hisada, Y., and S. Yamamoto. One-, two-, and three-dimensional site effects in sediment-filled basins.
Proceedings of the 11th world conference on earthquake engineering. 1996.
[2] T. Ichimura, M. Hori, J. Bielak. A hybrid multiresolution meshing technique for finite element threedimensional earthquake ground motion modelling in basins including topography. Geophysical JournalInternational, Volume 177, Issue 3, Pages 1221 1232, June 2009.
[3] Ichimura, T., Fujita, K., Errol, P., Quinay, B., Maddegedara, L., Hori, M., Tanaka, S., Shizawa, Y.,
Kobayashi, H., & Minami, K. (2015). Implicit Nonlinear Wave Simulation with 1.08T DOF and 0.270T
Unstructured Finite Elements to Enhance Comprehensive Earthquake Simulation. SC'15. IEEE, 2015.

36

[4] 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[5] Melosh, H., & Raefsky, A. (1981). A simple and efficient method for introducing faults into finite
element computations. Bulletin of the Seismological Society of America, 71(5), 1391–1400.
https://pubs.geoscienceworld.org/ssa/bssa/article-abstract/71/5/1391/118209/A-simple-and-efficientmethod-for-introducing?redirectedFrom=fulltext[6] Lysmer, J., & Kuhlemeyer, R. L. (1969). Finite dynamic model for infinite media. Journal of theEngineering Mechanics Division, 95(4), 859-878.
[7] Cerjan, C., Kosloff, D., Kosloff, R., & Reshef, M. (1985). A nonreflecting boundary condition for
discrete acoustic and elastic wave equations . Geophysics, 50(4), 705–708.
http://doi.org/10.1190/1.1441945[8] 林宏一, 引間和人. 差分法による三次元粘弾性波動場計算, 物理探査学会学術講演会講演論文集, 25-29,
103, 2000.
[9] Hisada, Yoshiaki. An efficient method for computing Green's functions for a layered half-space withsources and receivers at close depths. Bulletin of the Seismological Society of America 84.5 (1994):
1456-1472.
[10] Anderson, J. G. (2004). Quantitative measure of the goodness-of-fit of synthetic seismograms.
Proceedings of the 13th World Conference on Earthquake Engineering, 243.
https ://www.iitk.ac.in/nicee/wcee/article/13_243.pdf[11] Sanchez-Sesma, Francisco J., and Francisco Luzon. Seismic response of three-dimensional alluvial
valleys for incident P, S, and Rayleigh waves. Bulletin of the Seismological Society of America 85.1
(1995): 269-284.
[12] 林田智宏, 竹中博士 , 岡元太郎: 速度-応力型スタガード格子差分法を用いた 2 次元及び 3 次元地震波動
計算コードの作成, 九州大学理学部研究報告 地球惑星科学, 20 巻, 3 号, pp. 99-110, 1999.
[13] Graves, R.: Simulating seismic wave propagation in 3D elastic media using staggered-grid finite
differences, Bull. Seism. Soc. Am., Vol. 86, No. 4, pp. 1091-1106, 1996.
[14] Clayton, R., Engquist, B.: Absorbing boundary conditions for acoustic and elastic wave equations,
Bull. Seism. Soc. Am., Vol. 67, No. 6, pp. 1529-1540, 1977.
[15] Zeng, C., Xia, J., Miller, R., Tsoflias, G.: An improved vacuum formulation for 2D finite-difference
modeling of Rayleigh waves including surface topography and internal discontinuities, Geophysics, Vol.
77, No. 1, pp. T1-T9, 2012.
[16] Pitarka A.: 3D Elastic finite-difference modeling of seismic motion using staggered grids withnonuniform spacing, Bull. Seismol. Soc. Am., Vol. 89, No. 1, pp. 54-68, 1999.
[17] 笠松健太郎, 山中浩明: 中小地震の地震動シミュレーションによる 2004 年新潟県中越地震震源域の深部
地盤構造モデルの構築, 物理探査, 59 巻, 5 号, pp. 475-484, 2006.
[18] 笠松健太郎, 加藤研一: 大規模並列計算による 2011 年東北地方太平洋沖地震(Mw9.0)の長周期地震動シ
ミュレーション, 地球惑星科学連合大会, 2016.
[19] 吉村智昭, 永野正行, 久田嘉章, 青井真, 早川崇, Citak, S., 松島信一, 大西良広: 強震動予測手法に関するベンチマークテスト：数値解析手法の場合(その 1), 日本建築学会技術報告集, 17 巻, 35 号, pp. 67-72,
2011.
37

[20] Ichimura, T., Hori, M., Quinay, P. E., Wijerathne, M. L. L., Suzuki, T & Noguchi, S. (2012).
Comprehensive numerical analysis of fault-structure systems – Computation of the large-scale seismic
structural response to a given earthquake scenario –. August 2011, 795–811.
https://doi.org/10.1002/eqe.1158[21] Okada, Y. (1992). Internal deformation due to shear and tensile faults in a half-space. Bulletin of the
Seismological Society of America, 82(2), 1018–1040. https://pubs.geoscienceworld.org/ssa/bssa/articleabstract/82/2/1018/119580/Internal-deformation-due-to-shear-and-tensile[22] Fujita, K., Ichimura, T., Koyama, K., Inoue, H., Hori, M., & Maddegedara, L. (2017). Fast andScalable Low-Order Implicit Unstructured Finite-Element Solver for Earth’s Crust Deformation
Problem. Proceedings of the Platform for Advanced Scientific Computing Conference, 11.
[23] Gharti, H. N., & Langer, L. (2019). Spectral-infinite-element simulations of coseismic and postearthquake deformation. 1364–1393. https://doi.org/10.1093/gji/ggy495[24] Koketsu, K., Miyake, H., & Tanaka, Y. (2009). A proposal for a standard procedure of modeling 3-D
velocity structures and its application to the Tokyo metropolitan area , Japan. Tectonophysics, 472(1–4), 290–300. https://doi.org/10.1016/j.tecto.2008.05.037[25] Koketsu, K., Miyake, H., & Suzuki, H. (2012). Japan integrated velocity structure model version 1.
Proceedings of the 15th World Conference on Earthquake Engineering, 1773.
[26] Agata R, Ichimura T, Hirahara K, Hyodo M, Hori T, Hashimoto C, Hori M (2015) Numericalverification criteria for coseismic and postseismic crustal deformation analysis with large-scale highfidelity model. Proc Comput Sci 51:1534–1544. https ://doi.org/10.1016/j.procs .2015.05.344[27] Hori, T., Agata, R., Ichimura, T., Fujita, K., Yamaguchi, T., & Iinuma, T. (2021). High-fidelity elastic
Green’s functions for subduction zone models consistent with the global standard geodetic reference
system. Earth, Planets and Space, 1–19. https://doi.org/10.1186/s40623-021-01370-y[28] 富田史章，飯沼卓史，縣亮一郎，堀高峰，市村強，藤田航平，山口拓真 (2020). Trans-dimension インバージョン手法による沈み込み帯における固着すべり分布推定, 日本地震学会 2020 年度秋季大会.
（２）地震に関する災害被害予測のための大規模アプリケーションの開発
本項目では、「富岳」の性能を引き出すように、計算科学・計算機科学の最先端技術を駆使して、地震に関する災害被害予測のための大規模アプリケーションを改良することを目指している。具体的には、地震に関する災害被害予測のためのアプリケーション群の基幹数値解析コンポーネントである、地震・地震動・地殻変動計算等を行う大規模数値解析用コード群に対し、ソルバーの前処理等の改良を行い、「富岳」
実機上でより良好な性能を発揮するように開発を行うとともに、大規模数値解析コード群により解析する、全球モデルを視野に入れた超大規模解析モデルの構築手法を解析手法の特性を考慮しつつ開発することとしている。
本年度は、より効率的な大規模数値解析を可能とするように「①「富岳」の計算機特性を考慮した「富岳」の性能を引き出すソルバーの開発」、構造モデルの詳細度にあわせて超大規模モデルのデータ構造を変更した効率的な解析手法として「②ニューラルネットワークによりデータ構造を組み換えた効率的な前処理の開発」を行った。また、開発している大規模数値解析コード群の妥当性・有効性を検討するために「③開発された大規模数値解析コード群の実展開を踏まえた解析手法の有効性の確認」を行った。以下に、上記の①～③の具体的な成果をまとめる。
① 「富岳」の計算機特性を考慮した「富岳」の性能を引き出すソルバーの開発
38

①-(a). はじめに地震はときとして致命的な被害を人類社会へ引き起こすため、その被害軽減を目的としたさまざまな検討が行われてきた。検討のうち代表的なアプローチの一つは物理法則に基づく数値シミュレーションによる現象の理解と対策の立案であり、プレートの運動による地殻変動、地震発生から地震動伝播、都市の応答といった地震の各物理過程に対してシミュレーション手法が提案されてきた。これらの解析対象は広域かつ構造物に被害を及ぼす周波数帯を勘案すると必要な分解能が非常に細かくなるため、膨大な解析コストを解決し詳細な応答の検討が可能となるよう数値シミュレーションの高度化が図られてきた。そのうちのひとつとして、構造物被害の軽減を目的とした、構造物の地震応答解析の高度化がある。昨今では空間情報の蓄積・IoT や 5G による計測データの増加に伴い、従来よりも格段に高詳細な数値解析モデルを構築することが可能となりつつあり、特に超重要構造物の地震応答解析でこのような数値解析モデルを用いた
地震応答解析の実現が望まれている。
このような地震応答解析では、102-3 x 102-3 x 102 m 程度の領域で 10-2 m 程度の分解能で定義された高詳細な地盤・構造物系を対象とする。Geometry・局所的な非線形現象を正確にモデル化すること、応力フリーの境界条件を満足させることが必須であることから低次非構造有限要素による非線形動的解析が適している。また、Geometry を正確にモデル化する際に局所的に生成される小さな要素による厳しすぎる
Courant 条件を回避するために陰的時間積分が適している。そのため、従来の構造物地震応答解析と同様に、三次元低次非構造有限要素による陰的非線形動的解析を用いることとなるが、上記のような高詳細なモデルを対象とした場合、その数値解析は超大規模なものとなる。例えば、本研究における適用例では四面体二次要素を用いた有限要素モデルを用いているが、その自由度は約 490 億となる。そのため、陰的非線形動的解析で現れる線形方程式を Ax=f とすれば（非線形解析の詳細は section①-(b)参照）、この求解を直接解法で行うことはスーパーコンピュータを用いたとしても現実的ではない。また、疎行列とはいえ A
をメモリ上に確保すること自体が危ぶまれる状況にあるため、A をメモリ上に陽に確保するのではなく、反復解法の使用を前提として行列ベクトル積 Ax を ΣAexe のように要素剛性マトリクスをオンザフライで構築する element by element (EBE)法が適していることとなる。すなわち、反復解法である前処理付き共役勾配法を用い、その中のメインカーネルである Ax を EBE 法で置き換えることが大規模問題を解く上で適した方法となる（この詳細は section①-(b)参照）。しかしながら、低次非構造有限要素を用いた EBE 法
は、このような地震応答解析だけでなく多くのアプリケーションで現れる「演算量に対するデータの読み
書き量が多い」カーネルの典型例となっており、昨今の計算機構と相性の悪い方法であり、昨今の計算機構の性能を引き出した解析速度を達成するためには計算機科学的な検討を踏まえた特別なケアが必要となる（このようなカーネルと昨今の計算機アーキテクチャに関する検討は section ①-(c)参照）。事実、本研
究では計算機科学的な検討を踏まえた場合とそうでない場合で、実行時間に 5-6 倍の開きが生じており、このような検討の有用さを示している。なお、当たり前の話ではあるが、解析手法の近似で高速化されているわけではなく、計算機の機構を踏まえた適切な実装方法を開発し実装しているだけなので、得られる解は同じである。解析の高速化は計算科学の主眼のひとつであるが、昨今の計算機構の先鋭化に伴い、本研究のような計算科学をさらに加速する計算機科学的検討の重要性はさらに増していくと思われ、本研究のような検討は多くのアプリでみられる「演算量に対するデータの読み書き量が多い」カーネルの高速化に示唆を与えるものと期待される。
以上を踏まえ、①「富岳」の計算機特性を考慮した「富岳」の性能を引き出すソルバーの開発を行った。以下に①の構成を示す。まず、本研究で用いる低次有限要素による陰的非線形動的解析手法の詳細を
section①-(b)にて示す。なお、本研究では、従来法として前処理付き共役勾配法（PCGEBE）を、超大規模解析を可能とする超並列が可能な state of art のひとつである解の疎密・精度混合演算を用いた共役勾配法(SC14solver)を検討の対象とするので、これらについても説明する。次に、section①-(c)にて、多くのア39

プリへの示唆が期待される「演算量に対するデータの読み書き量が多い」カーネルの典型例である EBE について、その特性を昨今の計算機構を踏まえて考察するとともに、昨今注目を浴びている最新のアーキテ
クチャである Arm SVE アーキテクチャを対象として実装法の検討を行う。Section①-(d)にて、PCGEBE及び SC14solver に対して、提案手法を実装しその性能を Arm SVE ベースの「富岳」上で詳細に分析することでその有効性を確認する。最後に section①-(e)において、本提案手法及びその有効性のまとめを示す
とともに、「演算量に対するデータの読み書き量が多い」カーネルに対する今後の計算機アーキテクチャの進歩の動向にあわせた展望を示す。
①-(b). ターゲット問題
1) ターゲット問題
本研究では、非線形物性を持つ非均質構造中の波動伝播を対象に以下の式を解く。
(2-1)
ここで、ρ, u, t, x, c, f はそれぞれ密度、変位、時間、座標、材料定数、外力である。 材料定数 c は位置に応じて変化し、さらに歪・応力履歴により時間方向にも変化する。本研究で対象とする地盤・構造物の連成問題は、複雑形状・非均質物性を持った問題となる。この問題の求解のため、本研究では空間方向に有限要素法・時間方向に Newmark-β 法を用いる。対象の方程式は以下となる。
(2-2)

(2-3)

ここで、δu, u, v, a, q, f はそれぞれ変位増分、変位、速度、加速度、内力、および、外力となる。M, C, Kは整合質量マトリクス、レイリー減衰マトリクス、及び、剛性マトリクスとなる。dt, n は時間増分と時間
ステップ番号となる。ここでは、複雑形状における応力・歪応答を低コストで求解できる四面体二次要素を用いる。非線形 3 次元動的変位応答 un の求解では以下を繰り返す：

40なお、step 2)においては任意の材料構成則を用いることが可能であるが、本研究においては地盤のモデル化においては修正 RO モデル・Masing 則を用いる。計算コストのほぼ全てが式(2-2)の求解となるため、この式の求解に適したソルバーを開発する。
2) 線形方程式ソルバーと EBE 法
式(2-2)の超大規模方程式の求解においては、スケーラビリティ・収束性・ロバスト性・メモリ使用量を加味してソルバーを選択することとなる。Section①-(a)で述べたように、メモリ使用量が少なくスケーラ
ビリティ・ロバスト性が高いスタンダートな従来法はブロックヤコビ前処理を使った共役勾配法に EBE 行
列ベクトル積を用いた PCGEBE となる。そこで本研究の有効性を示すためこの PCGEBE の EBE に本手法を適用する。また、PCGEBE をベースに効率化を図った state of art solver はいくつかあるが、今回はSC14solver [1]を対象として、本手法の有効性を示すこととする。なお、SC14solver から派生した SC15[2]、IPDPS [3]、SC18 [4] などがあるが、これらのアルゴリズムは SC14solver よりさらに複雑になっており、提案手法の有効性の真水の部分を読者に対してわかりやすく説明することが困難であるため、
SC14solver を題材としている。以下、PCGEBE と SC14solver のアルゴリズム、及び、これらのソルバーのコアカーネルとなる EBE について説明する。
3) PCGEBE ソルバー
低次非構造格子有限要素法で離散化した式(2-2)のマトリクス A は超大規模な疎行列となるため、式(2-2)
を直接法にて求解することはメモリ使用量上困難で、反復法を使うこととなる。マトリクス A は正定値対称行列となるため、一般的に共役勾配法系の反復ソルバーを使うこととなる。共役勾配法の前処理には多数の方法が提案されているが、超並列環境におけるロードバランス・スケーラビリティを勘案して、節点毎にローカルな逆行列を前処理行列とするブロックヤコビ前処理が広く使われている。式(2-2)は各節点において x,y,z の三方向の変位が未知数となるため、3 x 3 のブロックヤコビ行列を使うこととなる。
Algorithm 1 に前処理付き共役勾配法のアルゴリズムを示す。ここでは行列ベクトル積(line 11)・ベクトル
内積(line 6,12)・ブロックヤコビ行列とベクトルの積(line 5)・AXPY (αx+y) (line 10,14,15)の 4 種類の計算が行われており、対象メッシュの要素数・節点数をコア間で均等に分割することで、いずれも計算負荷が均一になる。有限要素法により得られる A の場合は、計算コストのほとんどが行列ベクトル積となるた
め、この計算で用いる EBE の高速化が重要となる。

41

4) EBE 法
近年の計算機においては、演算性能に対するメモリ容量およびメモリ転送性能が相対的に低下する傾向となっている。従って、メモリ使用量・メモリ転送量を減らすアルゴリズム開発が求解可能な問題規模の拡大、及び、実行時間の短縮に重要となっている。そこで、マトリクスフリーの行列ベクトル積である
EBE 法を行列ベクトル積 f ← A u の計算に用いる。ここで、u, f は変位及び節点力ベクトルであり、A は
要素行列 Aei を重ね合わせて生成した全体行列である。Aei は各節点の座標(x)及び要素物性から計算可能である。通常の行列ベクトル積においては f ← A u の計算においてはメモリ上に格納した全体行列 A をメモ
リから読みだして u に掛け合わせる手順をとる。その一方で、EBE 法においては要素毎にローカルな要素
左辺ベクトル(2-4)
を計算したのち、これを全体左辺ベクトルに(2-5)のように足しこむことで計算を実施する。ここで、Qei は要素節点番号を全体節点番号にマッピングする行列である。座標・変位・節点力ベクトル（x, u, f）のサイズは小さく、節点・要素のオーダリングにより座標・変位・節点力ベクトルにおいて頻繁にアクセスする部分をキャッシュ上に載せておくことが可能なため、全体行列をメモリに格納する通常の行列ベクトル積計算手法と比べて、メモリからの読み出し量を大幅に削減可能である。その一方で、EBE 法では式(2-5)において変数のランダム足しこみが生じるため、
メニーコア・SIMD 計算機における並列計算は straightforward ではない。従って、ターゲットアーキテ
クチャに即した EBE 計算手法を選択することが高速な PCGEBE を実現するために重要となる。
5) SC14solver
42

PCGEBE をベースに解空間と演算精度空間の粗密を利用して効率化を図った手法の一つが SC14solver
である [1]。このソルバーにおいては、可変前処理付き共役勾配法、マルチグリッド法、及び、精度混合演算が用いられている。共役勾配法の本体部分は倍精度、前処理部分においては単精度が用いられており、また、この前処理においてはマルチグリッド法により計算コストの低減が図られているものの、全ての構成要素は PCGEBE をベースに作られているため、その中で用いられる EBE 法の性能が SC14solver の性
能を引き出すために重要となる。
以上をまとめると、PCGEBE も SC14solver も省メモリでスケーラブルに設計されており、性能は EBE
行列ベクトル積部分によることとなる。しかしながら、この部分はランダムアクセスを多く含んだ計算となるため、並列計算機機構に応じて最適な計算アルゴリズム・実装が変化することとなる。本研究ではこの部分に着目して研究開発を実施する。
①-(c). Arm SVE ベースのシステム用の EBE カーネルの開発
本研究では、Section①-(b)において説明された EBE カーネルを例に、「演算量に対するデータの読み書き量が多い」カーネルの高速化手法を開発する。その特性を昨今の計算機構を踏まえて考察するとともに、昨今注目を浴びている最新のアーキテクチャである Arm を対象として実装法の検討を行う。
1) Arm SVE アーキテクチャと「富岳」
Arm アーキテクチャのなかでも本研究では、低消費電力・高パフォーマンス・低プログラミングコストの特性から今後普及が期待される、wide SIMD 型の CPU の一種である Arm scalable vector extension(SVE)アーキテクチャに着目する。Arm SVE ISA では、従来からの軽量な 64bit Arm 命令セットに、128
bit から 2,048 bit までの可変長の SIMD 命令を追加することで、低消費電力かつ高い浮動小数点演算性能を狙ったものとなる。実際、Arm SVE の最初の実装である富士通の 512-bit SIMD メニーコア A64FXCPU ベースのシステムは 2019/11 の Green 500 list において 16.8 GFLOPS/W で 1 位となったように低消費電力を実現し、A64FX CPU を搭載した理化学研究所の「富岳」は 14.6 GFLOPS/W という大規模システムとしてはかなり高い電力性能を実現している。この電力性能の高さは ISA や製造プロセスだけでなくマイクロアーキテクチャにもよるものであり、Xeon などの他のメニーコア・wide SIMD CPU とは異なった特性を持つ。以下、その特性に応じたアルゴリズム・実装開発を説明する。
2) Arm SVE ベースのシステムにおける EBE カーネル
データアクセスの方法には計算機毎にいくつかの方法があり、また、それらのコストも計算機毎に異なることとなる。そこで、CPU と GPU など計算機機構に応じて異なる足しこみ手法がとられてきた。このように、この部分のアルゴリズムはハードウェアに依存して変更する必要が生じる。図(2).1 にマルチコア
CPU 用の標準的な EBE 計算アルゴリズムを示す。ここでは各コアにてテンポラリ配列 ftmp を確保し、lines 5-10 にて初期化したのち、各コアにてそれぞれの結果が lines 11-34 において ftmp に足し込まれる。最後に、コア毎の結果が全体左辺ベクトル f に足しこまれる（lines 46-56）. この方法は SC14solverにて用いられ、マルチコア CPU ベースの「京」コンピュータにおいて高い性能を発揮した。しかしながら、lines 31-33 におけるランダムデータ足しこみが最内となる ie ループの SIMD 並列化を阻害しており、wide SIMD CPU においては性能の著しい低下が予想される。さらに、メニーコア CPU においてはftmp の使用メモリ量も大きくなることとなる。そこで、近年の CPU において高性能を発揮するためには
EBE カーネルの計算アルゴリズム開発が必要となる。

43

図(2).1 As is EBE kernel algorithm for the K computer. Here we show the case for linear
tetrahedral elements with 4 nodes, each with 3 degrees-of-freedom corresponding to the x, y, and z
directions.
本研究では、まず SIMD を有効活用するアルゴリズムを開発した(図(2).2 参照)。ftmp へのランダム足し
こみが SIMD 演算器の適用を阻害していたため、ここでは最内ループを fie ← Aie uie の計算部(図(2).2 lines
17-32)、及び、要素左辺ベクトル fie の全体左辺ベクトル f への足しこみ部 (lines 34-42)の二つのループに分割している。これにより、計算量の多い前半部に SIMD 計算を適用できるようになる。ここでは、大きさ NL のループブロッキング(ここで NL は通常 SIMD 幅と同じ長さに設定する)を用いることで Aeu11-34などのテンポラリバッファをキャッシュ上にとどめることが可能となる。

44

図(2).2 Developed EBE kernel algorithm with SIMD enhanced computation and thread
partitioning
さらに本研究では、各スレッドに割り当てられるテンポラリ配列 ftmp を廃止するために、スレッド分割方法を開発した。通常のカラーリング方法では全メッシュにおいては各カラーで節点を共有しない要素が選ばれるのに対し（図(2).3）、提案手法においてはあらかじめ定めたスレッド数に応じて要素分割を実施す
る(図(2).4)。ここでは、グラフ分割手法(METIS)を用いて全体メッシュをスレッド数で分割し、節点を共有する要素を取り除くことで要素分割を実施している(図(2).4 の Color #1)。残った要素は再帰的に分割することで全要素を分割する(図(2).4 の Color #2, #3)。提案手法は節点データ u, x, f の再利用率が上がり、高い性能につながると期待される。

45

図(2).3 Standard coloring method. Here, elements in a single color are computed in parallel among
threads, and threads are synchronized among colors. While temporary buffers can be avoided,
nodal values (u, f) cannot be reused on cache.

図(2).4 Developed thread partitioning method. By decomposing the mesh with graph partitioning
methods, nodal values (u, f) can be reused more efficiently when compared to the standard coloring
method.
①-(d). 数値実験
この section においては、PCGEBE 及び SC14solver に対して、提案手法を実装しその性能を詳細に分析することでその有効性を確認する。ここでは、超詳細超大規模構造物モデルと収束特性が類似する weakscaling モデルセットとして、複雑地形を持つ地盤構造における非線形波動伝播問題を例に性能を計測する。性能計測は Arm SVE CPU スパコンである「富岳」と前世代の CPU スパコンの一つとなる「京」コンピュータの二種類のシステム上で実施する。
1) 問題設定
本研究では、図(2).5 に示すような 2 層の地盤構造を模擬した有限要素モデルを用いて非線形波動場解析に係る時間を計測した。第一層は非線形物性を持つ軟弱層であり、2 層は線形な基盤層である。このモデル
を x, y 方向に複製することで、表 1 に示す weak scaling 用の計測モデルを構築している。これらのモデルは全て最小要素サイズ 0.44 m という高分解能で離散化されており、非線形化した地盤においても数値解の収束を見込むことができる。ここでは、神戸海洋気象台において 1995 年の兵庫県南部地震時に観測された
神戸波 3 時間ステップ(dt = 0.002 s)の求解に要した時間を計測した。

46

図(2).5 Performance measurement setting. Minimum element size of 0.44 m is used.
表(2).1 Configurations of the models used to measure performance. A-Fugaku is a full scale model
for full Fugaku, while A-K is the same problem on the full K computer.

2) 計算機環境
本研究では、「富岳」、及び、「京」コンピュータにおいて性能を計測・比較した。
「富岳」は 48 計算コアを有する Arm v8.2-A SVE ベースの A64FX CPU 1 個を搭載する計算ノード
158,976 台からなるシステムである。各コアの SIMD 演算器は 512 ビットであり、各ノードのピーク性能
は倍精度で 32 × 48 コア × 2 GHz = 3,072 FP64 GFLOPS となる(2.2 GHz のクロックブーストで 3,379
FP64 GFLOPS)。FP32 は FP64 のハードウェアピークの倍となる。各計算ノードは HBM2 メモリ 32 GBを搭載し、ピークメモリバンド幅は 1,024 GB/s となる。A64FX CPU は 12 コア(core memory group;
CMG)毎にまとまった NUMA アーキテクチャとなるため、本研究では、4 つの CMG においてそれぞれ
MPI プロセスを 1 個立ち上げ、そこからそれぞれ 12 OpenMP スレッドを立ち上げる方法で並列計算をした。また、クロック周波数は 2.2 GHz で計測を実施している。
「京」コンピュータは 8 コアを有する SPARC64 VIIIfx CPU 1 個を搭載する計算ノード 82,944 台からなるシステムである。各コアの SIMD 演算器は 128 ビットであり、各ノードのピーク性能は倍精度で 8 ×
8 コア × 2 GHz = 128 FP64 GFLOPS となる。FP32 は FP64 のハードウェアピークと同じである。各計算ノードは DDR3-SDRAM メモリ 16 GB を搭載し、ピークメモリバンド幅は 64 GB/s となる。本研究では、各 CPU において MPI プロセスを 1 個立ち上げ、そこから 8 OpenMP スレッドを立ち上げる方法で並列計算をした。
3) EBE カーネルの性能
まず、開発した EBE カーネルの性能を「富岳」で計測した。ここでは、四面体二次要素カーネル(FP64)の性能を計測する。As is は図(2).1 における実装であり、[1]において用いられているものである。図(2).647

から、この as is カーネルの性能は必ずしも高いわけではないことがわかる(FP64 ピーク比 2.63%)。これは SIMD 命令の比率が 16.3%と低いことが一因であり、提案手法においてこの比率を 58.3%まで高めることで 3.06 倍という高速化につながっていることがわかる。さらに、スレッド分割手法によりメモリアクセ
ス量を削減することで、14.1 s (as is の 3.68 倍速)に、さらに、ループ分割による spill fill の削減と命令ス
ケジューリングの改善により 7.9 s になることがわかる。結果として、提案手法は as is の 6.56 倍速(FP64
ピーク比 13.1%)という高い性能が得られていることがわかる。

図(2).6 EBE kernel performance on one node of Fugaku. The measured elapsed time for conducting
one matrix-vector product in FP64. Numbers in brackets indicate ratio of attained FP64 FLOPS to
hardware peak.
4) PCGEBE における提案手法の性能
EBE のカーネル性能がほぼそのままソルバーの性能に直結する PCGEBE ソルバーにおける本手法の効果を「富岳」上で測る。カーネル単体での性能計測時の傾向と同じく、各アルゴリズム開発により高速化が達成され、最終的に as is 比で 5.1 倍速になっていることがわかる。この際の実行効率は 12.7%となり、低次非構造の陰解法有限要素法としては高い値になっていることがわかる。

図(2).7 Comparison of methods on model A-1 on 288 nodes of Fugaku. The measured elapsed time
for solving three time steps. Numbers in brackets indicate ratio of attained FP64 FLOPS to
hardware peak.

485) Effectiveness of developed method on SC14solver
続いて、state of the art solver の一例である SC14solver 中の EBE に本手法を実装し、性能を測る。
SC14solver では、前処理における単精度の EBE 計算が主要なカーネルとなる。「富岳」においては SIMD
幅が倍精度 8→単精度 16 に増えるので、PCGEBE における倍精度 EBE の場合に比べてハードウェアのSIMD 幅の増加分だけ提案手法における SIMD 計算の効果が見えるようになると期待される。実際、図(2).7 のように、PCGEBE では高速化率が 5.1 倍だったのに対し、SC14solver においては as is の 35.9 s(1.64% of peak)から提案手法は 4.13 s (13.5% of peak)と 8.69 倍速高速化していることがわかる。ArmSVE ISA は 2,048 bit まで対応しているなど今後も SIMD 幅が増えると予想されており、このような EBE
計算手法がさらに性能を出すために重要になると想定される。
本研究で対象とする重要構造物の詳細解析は大規模問題となるため、より多数ノードを用いて計算することとなる。その際、ノード数と問題規模を増やして計算する際の効率を示す weak scaling が重要となる。そこで、MPI プロセス当たりの問題規模を 164 万自由度に固定しつつ計算ノード数を増やすことでweak scaling を計測した。図(2).8 に計測結果を示す。 最小のモデル A-1 は「富岳」288 ノード上で 4.13
s (132 TFLOPS; FP64 ピーク比 13.5%)、最大のモデル A-Fugaku は「富岳」ほぼ全系に相当する 147,456
ノード上で 5.54 s (49.7 PFLOPS; FP64 ピーク比 9.97%)で求解され、288 ノードから 147,456 ノードへの
weak scaling efficiency は 74.4%と高い性能となった。

図(2).8 Sizeup scalability for model set A. The measured elapsed time for solving three time steps
using SC14solver. Blue: solver with as is EBE kernel, Orange: solver with developed EBE kernel.
6) 「京」コンピュータと「富岳」における性能比較
最後に、得られた「富岳」上での計算性能を「京」における性能と比較することで、計算機開発・計算手法開発による高速化率の内訳を確認する。本研究で用いた SC14solver の全系での性能計測問題(AFugaku)と同じものを「京」全系(A-K)で実行した際の time-to-solution と実行性能を表(2).2 に示す。
「京」全系での as is コードの実行時間は 328.0 s (FP64 ピーク比 8.60%)だったのに対し、提案手法は「富岳」ほぼ全系上で 38.8 s であった。
「富岳」と「京」のハードウェアピーク性能比は 46.9 倍であるため、
as is コードをそのまま使うだけではハードウェアの性能向上の恩恵の一部しか受けられないことがわかる。本研究の手法開発により、対「京」コンピュータ全系での性能倍率は「富岳」ほぼ全系において 59.2
倍となっている。

49

表(2).2 Performance of SC14solver on full K computer (82,944 node; peak performance 10.6
PFLOPS) and nearly full Fugaku (147,456 nodes; peak performance 498.3 PFLOPS [46.9-fold ofK])

このように、開発手法は標準手法である PCGEBE 及び発展的なソルバーである SC14solver の双方において「富岳」上で 5〜8 倍の高速化を実現しているなど、有効であることがわかる。また、「京」での as isの実行性能と比較することで、本研究のような手法開発が非構造格子有限要素法などの複雑なアプリケー
ションにおいてアーキテクチャの先鋭化に対応しシステムの性能を引き出すために必要となっていることがわかる。
①-(e). おわりに本研究では、重要構造物の超高詳細・超大規模有限要素地震応答解析を実施する際に必要となる、マト
リクスフリーの行列ベクトル積計算方法である EBE 法の高速化手法を提案した。ここでは、低消費電力な
Arm アーキテクチャのうち、高性能計算に着目した拡張がなされた Arm SVE ISA を用いて開発された「富岳」用の EBE カーネルアルゴリズムを開発した。
「富岳」での計測において、開発したカーネルを一
般的な前処理付き共役勾配法ソルバーに適用した場合で 5.11 倍速、大規模有限要素解析用ソルバーのstate of the art の一つである SC14solver に適用した場合で 8.69 倍の高速化を達成した。開発手法を用いることで、例えば、図(2).9 に示すような石油精製プラントなどの重要構造物の超詳細・超大規模解析も可能である。この場合での性能は、前世代の CPU スパコンである「京」上で実行された先端手法に対して、12.3 倍という大幅な高速化を達成しており、「富岳」の高い潜在性能と本研究の高性能カーネルプログラムを合わせることで 490 億自由度規模の超複雑プラントモデルにおいて 60 s 程度の実用的な時間スケールの地震応答解析も半日程度の時間で実行できるようになるなど、computational science の一つの主題である計算の高速化において計算機アーキテクチャに即した計算機科学検討を加えることで従来手法では解析することができなかった問題も解けるようになったことがわかる。
本研究で対象とした低次非構造有限要素における EBE は、昨今の計算機構と相性の悪い「演算量に対するデータの読み書き量が多い」カーネルの典型例ととらえられ、地震応答解析だけでなく多くのアプリケーションで現れる他の類似のカーネルにおいても有効な知見となると期待される。Arm SVE は仕様上
2,048 bit までの超 wide SIMD まで対応しているように、今後も一命令あたりの実行演算数を増やすこと
でハードウェアのピーク性能・電力性能の向上が進むと想定される。このようなアーキテクチャにおいては本研究のような計算機構の性能を引き出す計算機科学的な検討を踏まえたケアを実施することで、昨今の先鋭化した計算機の性能を十分に引き出しアプリケーションの性能向上と解ける問題の拡大につなげられると期待される。

50

図(2).9 Application problem results. Computation example of 49,064,764,344 degrees of freedom;
11,321,249,889 element model of fully coupled ground and structure earthquake shaking analysis
using developed solver on 20,736 nodes of Fugaku: (a) The 262 × 262 m domain was modeled with
1.6-cm tetrahedral elements with three soil layers and complex structure including pipes,
staircases and handrails (b,c). (d) The complex nonlinear seismic response was computed as aresult of the three-dimensional ground and structure configuration.

51

② ニューラルネットワークによりデータ構造を組み換えた効率的な前処理の開発②-(a). はじめに
非線形時間発展問題の応答評価の重要性は分野を問わず非常に高く、多くの分野でこの実現・効率化へ向けた研究がなされている。従来から超大規模解析の state of the art が集う Supercomputing Conferenceの Gordon Bell Prize Session においても partial differential equation (PDE)を解く非線形時間発展問題求解はメインテーマの一つであり、スパコン全系による複雑な形状を精緻にモデル化した 1 兆自由度を超える implicit solver による超高分解能非線形時間発展問題の応答評価なども実現されている(例えば[2])。その一方で、近年では現代の計算機アーキテクチャにおいて効率的に計算可能なニューラルネットワーク（Neural Network:NN）をはじめとする Machine Learning 手法が発展してきており、従来の PDE 求解のスキームに NN を加えることによる更なる解析能力向上を目指した研究もなされつつある(例えば[4])。
このような多様な研究がなされているものの、「解析コストは高いものの高精度で解を求めることを指向す
る従来の PDE 求解スキーム」と「解析コストを抑制できるが精度がやや低くなる傾向にある NN ベースの評価」は、そもそもの両者の出自・目的が異なることによる解の要求精度のずれのため、融合不足が生じている状況にあり、さらなる両者融合の真価を発揮する余地がある。そこで本研究においては、PDE の一つである連続体の波動方程式の求解スキームにおいて NN を活用できるように前処理・問題構造のデータ構造を変更することで、解析の効率化を図る手法を開発した。ここでは、PDE の特性を反映するグリーン関数を近似する NN (GF-based NN)を前処理で用いることで、陰解法ソルバーを高速化する方法を提案する。GF-based NN によりソルバーの収束性が改善するだけでなく、PDE の求解で生じるランダムメモリ
アクセスベースの計算を NN における連続メモリアクセスベースの高密度な計算に置き換えることで演算効率も向上し、従来のソルバー比で 4.26 倍の高速化を実現した。また、前処理計算はローカルに実施できるため超並列環境においても高い性能となり、「富岳」の 864 ノードから 55,296 ノードまでの weakscaling で効率 98.1%を達成した。以下、section ②-(b)にて提案手法の概要を示し、section②-(c)にて数値実験を通して提案手法の詳細を説明する。最後に、section②-(d)にて本研究のまとめを示す。
②-(b). 提案手法
PDE ベースの非線形動的解析においては、各時間ステップにおいて対象の PDE は線形化・離散化した線形方程式Aδu = fを解く。ここで、Aとfはnステップ目における既知のマトリクスとベクトルであり、これらの自由度が大きいため、反復解法を使ってδuを求解することとなる。本研究では、グリーン関数を模擬する NN を反復解法の前処理に用いることで、ソルバーの収束性改善を図る。以下、本 section において
は提案手法の概要を示し、section②-(c)において実際の問題における数値実験を通して詳細を説明する。
本研究では、A−1 を近似するBを従来の PDE 求解手法に比べて低解析コスト・低適用コストで実現する方法を開発し、前処理行列として用いる。Bの高次モードにおける近似性能がソルバーの収束性改善のために重要となるが、高次モード、及び、低次の transit なモードの双方を高精度で含む NN を低コストで構築することは簡単ではない。例えば、(x, Ax)の組のデータを多数保存し、このデータを活用して構築したNN を前処理に用いることが考えられるが、実際には保存できるデータの組の数は限られるため、コストの割に性能は上がらないこととなる。そこで本研究では、保存可能なデータ量で系の特性を高次モードまで捉えられるよう、対象の PDE の特性を反映するグリーン関数を介した NN を構築した。ここでは、データから直接 PDE の応答を学習するのではなく、データからグリーン関数を介して PDE の応答を予測する。一般に、PDE: 𝐿𝐿(𝑎𝑎(𝑥𝑥)) = 𝑏𝑏(𝑥𝑥) (ここで𝐿𝐿(𝑥𝑥)は linear differential operator、𝑏𝑏(𝑥𝑥)は既知関数、𝑎𝑎(𝑥𝑥): は系の応答)としたとき、デルタ関数に対する応答𝐿𝐿(𝑔𝑔(𝑥𝑥, 𝑠𝑠)) = 𝛿𝛿(𝑠𝑠 − 𝑥𝑥)をグリーン関数と呼び、この関数を用いることで𝑎𝑎(𝑥𝑥)は明示的に

52

(2-6)𝑎𝑎(𝑥𝑥) = ∫ 𝑔𝑔(𝑥𝑥, 𝑠𝑠)𝑏𝑏(𝑠𝑠)𝑑𝑑𝑑𝑑と求めることができる。マトリクス方程式Aδu = fとの比較からわかるように、これはδu = A−1 fの求解と等価であり、グリーン関数を求めることができれば PDE の求解に対応するA−1を求めることが可能である。

ところが、実際には𝑔𝑔(𝑥𝑥, 𝑠𝑠)は PDE の特性を反映して非常に複雑な分布を有し、非均質問題においては場所に応じて大きく変動することとなる。従って、𝑔𝑔(𝑥𝑥, 𝑠𝑠)を求める計算はAδu = fの求解と同等以上のコストがかかることとなり、これまで線形方程式の求解にグリーン関数はあまり用いられてこなかった。

本研究では、グリーン関数を用いた PDE の求解アプローチにおいて NN を使うことで、A−1を近似する

Bを低コストで求める。PDE の特性を反映するグリーン関数の分布は複雑となるが、一度高精度でグリーン関数を求めることができれば、PDE を高精度で求めることが可能となる。さらに、グリーン関数を少数のパラメータで表すことができるようなパラメタリゼーションにより、少ないデータから効率的に高精度な NN を構築できるようにした。次の section においては 3 次元波動方程式を対象に提案手法の詳細を説明する。
②-(c). 数値実験
GF-based NN は様々な反復解法の前処理に適用可能であるものの、本研究では共役勾配法への適用を通して手法の有効性を示す。非均質性の強い問題を対象に、計算機 1 台上での性能を計測したのち、多数計
算ノード上でのスケーラビリティを計測する。
1) 対象問題
本研究では、3 次元非均質物性領域における波動方程式を対象とする。この PDE の求解には各種の離散化方法を使うことができるが、今回は簡単のため空間方向に構造格子有限要素法（一次のボクセル要素）
を使い、時間方向には解の安定性が担保される Newmark-β法(β=1/4, γ=1/2)を使う。対象方程式はとなる。ここでA ⇐ (4𝑑𝑑𝑑𝑑 2

M+

2

𝑑𝑑𝑑𝑑

(2-7)

Aδu = fC 𝑛𝑛 + K 𝑛𝑛 )であり、f ⇐ b𝑛𝑛 − q𝑛𝑛−1 + C 𝑛𝑛 v 𝑛𝑛−1 + M(a𝑛𝑛−1 +

4

𝑑𝑑𝑑𝑑

v 𝑛𝑛−1 )である。ここで、δu, u, v, a, bはそれぞれ変位増分、変位、速度、加速度、外力であり、M, C 𝑛𝑛 , K 𝑛𝑛 はそれぞれ整合質量マト

リクス、レイリー減衰マトリクス、剛性マトリクスである。𝑛𝑛は時間ステップ番号、𝑑𝑑𝑑𝑑は時間増分である。
まず、式(2-7)を解き、その後以下のように変数を更新することで非線形動的解析を実施する。
q𝑛𝑛 ⇐ q𝑛𝑛−1 + K 𝑛𝑛 𝛿𝛿u𝑛𝑛

(2-8)

2
𝛿𝛿u𝑛𝑛
𝑑𝑑𝑑𝑑

(2-10)(2-9)

u𝑛𝑛 = u𝑛𝑛−1 + 𝛿𝛿u𝑛𝑛

v 𝑛𝑛 = −v 𝑛𝑛−1 +a𝑛𝑛 = −a𝑛𝑛−1 −4 𝑛𝑛−14v+ 2 𝛿𝛿u𝑛𝑛
𝑑𝑑𝑑𝑑
𝑑𝑑𝑡𝑡

(2-11)

ほとんどの計算コストは正定値対称な疎行列となるAに対する式(2-7)の求解にかかることとなるため、この高速化が重要となる。
上記の 3 次元非均質物性領域における波動方程式問題の例として、図(2).10 に示す頭部の CT スキャンを用いて作成した非均質モデルを対象に提案手法の性能を計測する。地震シミュレーションでは、インピー
53

ダンスコントラスト比が高い不均質な弾性体の動的問題を解くことになるため、このような性質を持つ地殻モデルによる開発手法のテストを当初計画した。しかしながら、論文の締切までに、論文に掲載可能な
このようなモデルの準備が難しいことが判明したため、等価な性質をもつ、公開データである頭部モデル
を用いて開発手法のテストを行うこととした。なお、複雑な不均質性の分布のみを頭部モデルから用いており、開発手法の肝となるグリーン関数の影響範囲を決定する物性は地震シミュレーションでよく用いられるものを用いてテストを行っているため、本開発手法の地震問題での有効性を検証するのに十分な問題設定となっている。この問題は Vs = 50-120 m/s の幅を持つ非均質な問題であり (Vp、密度 ρ、減衰率 hは全域で Vp=200 m/s、ρ=1000 kg/m3、h=0.1%で固定)、512 × 512 × 512 要素で離散化した際の自由度は
405,017,091 となる。時間増分には dt=0.001 s を用いる。weak scaling の計測のため、このモデルを x,y方向に複製し接続することで収束特性の似た問題セットを作成した。

図(2).10 Distribution of the Young's modules of the target head model2) GF-Based NN の構築式(2-7)の場合、グリーン関数 G が有意な大きさを取る範囲は対象領域の波速 Vp、Vs と時間刻み幅 dt によって決まり、3 次元の構造格子有限要素法において数値解の収束が保証されるような問題設定の場合は、対象地点を中心とした Nx×Ny×Nz 個の要素内に主要成分が含まれることになる。そのため、この限られた領域における入力 flocal と、その中心における解 ulocal を結びつけるグリーン関数 Blocal を NN として構築す
れば、精度の高い推定解を得られると期待される。本適用例の対象物性範囲では、Nx=Ny=Nz=6 とすれば良いことがわかったので、この条件の下で高次モードの発現をよくするために、式(2-7)の右辺にランダム項をいれて、Vs = 50-120 m/s の範囲内の様々な物性パターン（9,301 パターン）で数値シミュレーションを行い、flocal とその中心における推定解 ulocal を 1,620 万組以上生成し、これに対して誤差最小化の最適化を行うことで数値的なグリーン関数を求めている。なお、ここでは直接一般の物性のグリーン関数の形状を求めることはせず、2 段階の学習を経ることで精度改善を達成している。図(2).11 に構築した NN による
グリーン関数の分布を示す。上段が数値的に求めたグリーン関数の分布であるのに対し、下段は NN により推定されたグリーン関数と正確なグリーン関数の値の差となる。ここから、グリーン関数が非常に複雑な分布を持つこと、また、その複雑な分布が高精度で推定されていることがわかる。

54

図(2).11 Distribution of the Green's functions. The top row indicates the target values of theGreen's function distribution (G), while the bottom row indicates the difference (G-GNN)
between the target distribution and the distribution obtained using the generated NNs(figures appear to be rendered blank as the error is very small). The developed NNs
estimated the complex distribution of the Green's functions with a high accuracy for all of
the components shown. Similar accuracy was obtained for heterogeneous material.

3) GF-Based NN を活用した前処理開発した GF-based NN を共役勾配法の前処理として活用する（Algorithm 1）。境界付近ではグリーン関数の形が媒質内とは異なる特性で複雑に変化するため、今回の問題においては提案手法では領域内部をGF-Based NN で、境界付近は 3 × 3 ブロックヤコビ前処理を用いた共役勾配法で解く。また、MPI 領域境界においても 3 × 3 ブロックヤコビ前処理を用いた共役勾配法を用いる。3 次元ボリュームにおける境界
部・MPI 領域境界の割合は非常に小さいため、全解析コストに占める境界部の求解コストの割合は小さいと期待される。GF-based NN および境界部の更新からなる前処理は全て単精度演算で、その他の計算は全て倍精度演算で計算する。

55

まず、図(2).12 において、反復解法における誤差の収束履歴を示す。ここでは、提案手法による解と、3× 3 ブロックヤコビ前処理を用いた共役勾配法(CGBJ)において相対誤差 10-8 まで解いた際の解との差を示している。均質な部分の誤差が先に減り、その後反復が進むにつれて非均質部分の誤差が減っていることがわかる。5 反復目において誤差はランダムな分布になっており、誤差レベルが CGBJ によるものと同等となっていることがわかる。比較のために可視化した CGBJ の収束履歴をみると全域において均一に誤差が低減していることがわかる。CGBJ においては相対誤差 10-8 まで達するのに 42 反復要していたのに対し、提案手法は 5 反復しかかかっておらず、GF-based NN により B を高精度で求められていることがわかる。

56

図 (2).12 Distribution of the difference (u-ucorrect) between the partially solved solution andaccurate solution solved up to (ε=10-8) on the head model. Error indicate the accuracy of theobtained solution (relative error ||Au-f||2 / ||f||2 ). The range of the color scale was different for eachfigure as the norm of u-ucorrect changes significantly in each refinement iteration. Elapsed timemeasured on one Oakbridge-CX node.
次に、提案手法と CGBJ の実行時間を比較する。ここでは、東京大学情報基盤センターの OakbridgeCX を 1 ノード利用し、OpenMP により 56 コアを活用したスレッド並列プログラムを使って解析を実施した際の実行時間を比較する。提案手法・CGBJ ともに AVX-512 intrinsics を用いて SIMD 演算を明示的に指示してコードを構築しており、同等のチューニングとなっている。図(2).13 に各ソルバーの収束履歴を示す。CGBJ における反復当たりの実行時間は 1.40 s (FP64 ピーク比 10.9%)であったのに対し、提案手法の反復当たり実行時間は 2.3 s (FP64 ピーク比 24.1%)となった。反復あたりの実行時間は増えているものの反復数が大幅に減っているため、ソルバー全体で 4.26 倍の高速化となっている。なお、NN-based
method と PDE 求解スキームを組み合わせて精度保証をした解を得る簡易な方法である、NN-based
method による推定解を CG 法の初期解として使う方法では、CGBJ の 42 反復より少ない 36 反復で解を得られているものの、本研究のように解の改良に NN を使う方法と比べて性能向上幅が小さいことがわかる。

57

図(2).13 Convergence speeds of solvers for the head model measured on one Oakbridge-CX
node. The convergence ratio (ratio of change in relative error per iteration) was greatly
improved by the use of the GF-based NNs in each refinement iteration in the proposed
method. The simple scheme using GF-based NNs for initial solution of PDE solver scheme(CGBJ) could not utilize this improvement in the convergence ratio.
最後に、提案手法のスケーラビリティを「富岳」上で計測する。ここでは、Oakbridge-CX の Xeon
CPU 用に開発したコードの AVX-512 intrinsics を Arm の ACLE intrinsics に変更することで「富岳」のA64FX CPU において Xeon 機と同等のチューニング特性を持つコードを用いた。
「富岳」の各計算ノードにおいて 4 MPI プロセス × 12 OpenMP スレッド立ち上げ、CPU 周波数 2.2 GHz で性能計測を実施した。図(2).14 に計測結果を示す。GF-based NN 部、境界領域の更新部、その他の計算ともにスケールし、「富岳」の約 1/3 に相当する 55,296 ノードまで 98.1%の weak scaling 性能が得られた。55,296 ノードにおけるソルバー全体の性能は 30.3 PFLOPS (FP64 ピーク比 16.2%)となった。図(2).15 に示すように提案
手法は strong scaling 性能も良好であり、24 ノードから 768 ノードまでの strong scaling 性能は 66.7%となった。

図(2).14 Weak scaling results of proposed solver on Fugaku. The head model isduplicated in the x and y directions to make large models. Each compute node solves
a domain with 100.6 million DOF with inter-process communication among
adjacent domains. Numbers in brackets indicate ratio to FP64 peak.

58

図(2).15 Strong scaling results of proposed solver on Fugaku using a 2.4 billion DOF
problem constructed by connecting 6 copies of the head model.
通常、NN を PDE の求解に用いる場合には精度保証がなされないが、本研究で開発した手法は反復解法の前処理に用いているため解の精度が担保される。仮に、本研究で構築した GF-based NN を前処理に用いずに直接時間発展問題の解の推定に用いた場合、解の精度保証がなされないために図(2).16 のように多数
タイムステップ後の結果が発散してしまう(GF-based NN without accuracy assurance)。その一方で、提案手法は多数タイムステップ計算時においても通常の反復法ソルバーと同様に精度が担保された結果が得られていることがわかる。

図(2).16 Distributions of u at different time steps on the head model. The solutionof the proposed method with accuracy assurance was consistent with that of the
PDE solver scheme (CGBJ), while that of the NN-based method without accuracy
assurance (corresponds to terminating solver after one iteration of Algorithm 1)
diverged after 200 time steps (energy conservation of the input force was notsatisfied due to the limited accuracy). A typical PDE simulation would require tensof thousands of time-steps and thus accuracy assurance is required. Responsecomputed for oscillating gravitational force.
②-(d). おわりに本研究では、前処理・問題構造のデータ構造を変更することで効率的な前処理を達成する方法として、
非線形時間発展問題への NN の適用を念頭に微分方程式による解空間の拘束を効率的に行う NN を示しこ59

れを用いて、NN ベースの評価を任意精度で効率的に改善する収束性が保証された手法を提案した。ここでは昨今の計算機アーキテクチャへの親和性が低い従来の PDE 求解スキームを計算機アーキテクチャへの親和性が高い連続メモリアクセス卓越型の NN ベースの評価に解析コストを移行させることで、より高い実効性能で演算を可能とした。その結果本手法は、従来の PDE 求解のスキームである 3 × 3 ブロックヤコ
ビ前処理を用いた共役勾配法に比べ 5.48 倍高速に解を得ることに成功し、また、ローカルな演算のみからなる前処理手法としたために「富岳」の超並列環境においても高いスケーラビリティを達成した。
本研究では NN により微分方程式の構造であるグリーン関数を分析し、inferencing コストも抑えた NNを構築することに成功したが、もう少しパラメータ数を増やすことで NN の精度を向上させることは無論可能であり、解析コストのバランスを考慮して問題設定に応じた設定を選んで解析を行うことが可能であることも本手法の特徴である。本アルゴリズムは対称行列系を対象としているが、容易に非対称行列系に適用可能であり、より一般の非線形時間発展問題への適用が可能となる。また、今回は deterministic かつ
ノイズのない問題を適用例として示しているが、モンテカルロシミュレーションを行うことで、このような deterministic でない問題についても適用可能であり、今後の応用が期待される。
③ 開発された大規模数値解析コード群の実展開を踏まえた解析手法の有効性の確認③-(a). はじめに開発された大規模数値解析コード群の実展開を踏まえた解析手法の有効性を確認するため、大規模スパ
コン用に構築した計算手法を GPU 計算機へ移植し、その有効性を確認した。なお、実展開を踏まえ、より精度よく液状化等の現象を解析可能な構成則として、2018 PEER blind prediction contest で優勝したチームが使用している解析プログラム FLIP が採用している構成則を実装している。液状化解析などの構成則の計算は複雑で計算量が多いが、各要素に対して行う処理はほとんど同一という GPU が得意とする種類の計算である。そのため、GPU を使うことで、高速化が期待できる。GPU 環境でのプログラム開発はCUDA と呼ばれる専門の並列プログラミング言語が使用されてきた。しかし、CUDA での実装は開発コストが高く、かつ複雑になりがちであるため、多数の地盤・材料構成則の解析に CUDA を使用するのは現実的ではない。近年では、OpenACC と呼ばれる並列プログラミング言語の開発が進み、低開発コストでか
つメンテナンス性の高い、GPU 用プログラム開発が可能となってきている。例えば、[5]では、OpenACC
を使って大規模地震応答解析手法を GPU 環境へ移行し、CUDA での実装と比べても十分な性能を発揮している。これらを踏まえて、本研究では、OpenACC によって GPU 用の手法を開発し、小規模な計算環境
で、高詳細な 3 次元液状化解析が実行できるようにすることをめざした。開発した手法を GPU の搭載されたスーパーコンピュータ ABCI で性能計測を行い、既往手法では CPU ベースのス－パーコンピュータ上で行っていた解析が、本研究の手法を使うことで、小規模計算機クラスタ程度の計算資源で実施できるようになることを示した。また、適用例として、2 つの液状化解析を行った。1 つは、格子状地盤改良による既設住宅の耐震性強化を模した解析、もう 1 つは、河口付近の地形を模したモデルの解析である。本研究の手法によって、小規模な計算環境で、大規模な 3 次元地震応答解析が可能になったことにより、このような対策の性能照査もより容易に実施できるようになると考えられる。本研究は、ポスト「京」重点課題、「富岳」成果創出加速プログラムを通して開発してきた地震応答解析手法を実務で使用される構成則を使った解析に展開したものであり、これらの解析手法の実問題に対する有効性を示すものである。
③-(b). 手法
1) 対象問題
本研究では、2 次非構造要素の有限要素法による動的 3 次元地震応答解析を対象とする。有限要素法を使用する理由は、有限要素法が複雑な地盤や構造物の形状を正確にモデル化でき、また、分解能を自由に60設定できるため特定の部分だけを高分解能にすることができるため、特に大規模解析では計算量を削減できるからである。本研究の対象問題は非排水条件下の液状化解析であり、支配方程式は、以下の式(2-12)である。ただし、この式は、有限要素法と Newmark-β 法によって離散化されている。

𝑨𝑨 =(2-12)

𝑨𝑨𝛿𝛿𝒖𝒖(𝒏𝒏) = 𝒃𝒃

4
2
𝑴𝑴 + 𝑪𝑪(𝑛𝑛−1) + 𝑲𝑲(𝑛𝑛−1)
𝑑𝑑𝑡𝑡 2𝑑𝑑𝑑𝑑

𝒃𝒃 = 𝒇𝒇(𝑛𝑛) − 𝒒𝒒(𝑛𝑛−1) + 𝑪𝑪(𝑛𝑛−1) 𝒗𝒗(𝑛𝑛−1) + 𝑴𝑴 �𝒂𝒂(𝑛𝑛−1) +(2-13)
4 (𝑛𝑛−1)𝒗𝒗
�
𝑑𝑑𝑑𝑑(2-14)

ここで、𝑴𝑴, 𝑪𝑪, 𝑲𝑲はそれぞれ、質量行列、レイリー減衰行列、剛性行列である。また、𝛿𝛿𝒖𝒖, 𝒇𝒇, 𝒒𝒒, 𝒗𝒗, 𝒂𝒂 はそれぞれ、変位増分、外力、内力、速度、変位であり、𝑑𝑑𝑑𝑑 は時間ステップ幅である。∗(𝑛𝑛) は、変数∗が、𝑛𝑛 ステ

ップ目の変数であることを表す。以下の手順により、液状化解析を行う。
1. 外力𝒇𝒇を読み込み、係数行列𝑨𝑨と右辺ベクトル𝒃𝒃を計算する。
2. 式(2-12)を解き、変位増分𝛿𝛿𝒖𝒖をアルゴリズム 1 によって得る。
3. 以下の式により、変位𝒖𝒖、速度𝒗𝒗、 加速度𝒂𝒂を更新する。
𝒖𝒖(𝑛𝑛) = 𝒖𝒖(𝑛𝑛−1) + 𝛿𝛿𝒖𝒖(𝑛𝑛)𝒗𝒗(𝑛𝑛) = −𝒗𝒗(𝑛𝑛−1) +𝒂𝒂 = −𝒂𝒂(𝑛𝑛−1) −2𝛿𝛿𝒖𝒖(𝑛𝑛)
𝑑𝑑𝑑𝑑

44
𝒗𝒗 + 2 𝛿𝛿𝒖𝒖(𝑛𝑛)
𝑑𝑑𝑑𝑑
𝑑𝑑𝑡𝑡

(2-15)(2-16)(2-17)

4. 以下の式により、剛性行列 K と内力 q を更新する。

𝑲𝑲(𝑛𝑛) = � � 𝑩𝑩T �𝑫𝑫s
𝑒𝑒

𝑽𝑽𝒆𝒆

𝑒𝑒

𝑽𝑽𝒆𝒆

𝒒𝒒(𝑛𝑛) = � � 𝑩𝑩T �𝝈𝝈′

(𝑛𝑛)

(𝑛𝑛)

+ 𝑫𝑫w � 𝑩𝑩𝑑𝑑𝑽𝑽+ 𝑝𝑝(𝑛𝑛) 𝒎𝒎� 𝑑𝑑𝑽𝑽

(2-18)(2-19)

ただし、𝑩𝑩は変位ひずみ変換行列、𝑫𝑫𝑠𝑠 , 𝑫𝑫w , 𝝈𝝈′ , 𝑝𝑝はそれぞれ、地盤、水の弾塑性行列、地盤の有効応力、過剰間隙水圧である。これらは、構成則によって計算される。また、𝒎𝒎 = {1 1 1 0 0 0}であり、∫𝑽𝑽 ∗ 𝑑𝑑𝑽𝑽 は要𝒆𝒆

素𝑒𝑒での体積積分を表す。

本研究では、[6]による 2 次元の構成則を 3 次元に拡張した構成則[7]が使用されている。これは、実務で使用されており、2018 PEER blind prediction contest で優勝したチームが使用している解析プログラム
FLIP がもとにしている構成則である。この構成則は、マルチスプリングモデルの一種であり、以下の式のように、多数の 1 次元(1D)ばねの重ね合わせとして 3 次元構成則を表現している。

𝑁𝑁s

𝑫𝑫s = 𝑲𝑲u 𝒎𝒎𝒎𝒎T + � 𝑤𝑤𝑖𝑖s𝑖𝑖=1

d𝐹𝐹𝑖𝑖 (𝛾𝛾𝑖𝑖 )𝒏𝒏𝑖𝑖 𝒏𝒏T𝑖𝑖
d𝛾𝛾𝑖𝑖
61

(2-20)

𝑁𝑁s′𝝈𝝈′ = 𝜎𝜎m
𝒎𝒎 + � 𝑤𝑤𝑖𝑖s 𝐹𝐹𝑖𝑖 (𝛾𝛾𝑖𝑖 )𝒏𝒏𝑖𝑖(2-21)

𝑖𝑖=1′, 𝐾𝐾u はそれぞれ平均有効応力、体積弾性係数であり、𝛾𝛾𝑖𝑖 , 𝑤𝑤𝑖𝑖s , 𝐹𝐹𝑖𝑖 (𝛾𝛾𝑖𝑖 ), 𝒏𝒏𝑖𝑖 はそれぞれ、1 次元ばね𝑖𝑖の
ただし、𝜎𝜎m

ひずみ、重み係数、ばね力、3 次元ひずみを 1 次元ばねのひずみに変換するベクトル、である。構成則の評価に使用される𝑁𝑁𝑠𝑠 は 1 次元ばねの本数である。解の離散化依存性を避けるために、ばねの本数𝑁𝑁𝑠𝑠 は十分に多い必要があり、本研究では前述の FLIP に合わせて、𝑁𝑁𝑠𝑠 = 300とした。また、液状化解析では、局所的な液状化により、一要素内で物性が急激に変化するので、それを計算するために一要素につき 5 点で、構成則の評価を行っている。一つのばねにつき、7 つの変数（22 byte: 5 つの 32 ビット変数と 2 つの 1 ビ

ット変数）を保持する必要があり、また、1 つのばねの計算に約 100 実数演算が必要である。つまり、1 要
素当たり 33 KB のメモリと 1 タイムステップあたり 150,000 実数演算が必要である。例えば、後述する
適用例 2 で使われるモデルの非線形要素の数は、15,798,787 であるので、構成則のばねの情報を格納するのに必要なメモリ使用量は約 500GB、1 時間ステップ当たりの実数演算は 2.4 兆回であり、膨大なメモリと計算量が必要である。
以上をまとめたものをアルゴリズム 1 に示す。詳細は、次のセクションで説明する。

62

アルゴリズム１ 開発手法

63

アルゴリズム 2 前処理方程式の求解アルゴリズム（共役勾配法）

2) OpenACC と低精度変数と使用した GPU による計算本研究は、実展開を目指し複雑な地盤構成則を使った大規模な 3 次元地震応答解析を小規模な計算機クラスタ上でも実用的な時間内に実施できるようにすることを目的としている。小規模な計算機環境で、大規模な解析を短時間で行うには、演算密度を高める必要がある。本研究では、そのような計算が可能なハードウェアとして、近年、科学計算においても使われるようになってきている GPU を用いることとした。
GPU は数千のコアを搭載しており、超大量の並列計算を行うことによって高速に計算を行うことができる。GPU を使った計算の概要を図(2).17(a)に示す。図(2).17(b)に示されるように GPU は CPU に接続して使われる。CPU から GPU に並列計算の実行を指示することで、GPU 上で計算が実行される。ただし、
GPU は CPU メモリとは別にそれ自身のメモリを持っており、GPU 上での計算に必要なデータは GPU メ
モリ上にある必要があるため、GPU 上での並列計算に先立って、必要なデータを CPU から GPU に送る必要がある。また、GPU 上での計算の結果を必要に応じて CPU メモリに戻す必要もある。これらの実装が必要であるため、GPU 用のアプリケーション開発には、CPU 上のアプリケーション開発にはない難しさがある。GPU 上で、本研究が対象とする、複雑な構成則を用いた非構造要素の有限要素法による地震応
答解析を実行するにあたって考えるべきことが 3 つある。1 つ目は、GPU 用のコードの開発とメンテナン
ス性である。実務で使うことを考えると、できるだけ低開発コストで、かつ、メンテナンスがしやすい必要がある。2 つ目は、GPU に適したアルゴリズムである。GPU を使うことよって高速化されやすい計算とあまり高速化されない計算があるため、支配方程式の求解法や構成則の計算において、GPU の性能が発揮できるようなアルゴリズムを採用する必要がある。そして、3 つ目が GPU、CPU 上でのデータ管理である。CPU-GPU 間のデータ転送は比較的遅いので、データ管理を適切に行い、CPU-GPU 間のデータ転送を最低限にする必要がある。

64

図(2).17 GPU を使った計算
GPU 用のコードの開発には、従来は、デファクトスタンダードである CUDA と呼ばれる専用の並列プ
ログラミング言語が使われていた。地震応答解析に関しても、CUDA を使って GPU 環境へ移行し解析を高速化するという研究も行われてきている（例えば[5][8]）。しかし、性能が高いコードを実装するには、GPU 上で動くコードを使って実装するだけでなく、GPU のアーキテクチャに関する知識も必要であり、CUDA を用いた場合の開発コストは膨大となる。また、CUDA を使って開発したプロクラムは CPU 上では動かず、特定の GPU 上でしか動作しない。さらに、開発されたコードは実務で継続的に使用していくものであり、将来のメンテナンスや修正、アップデートが難しくなってしまう。本研究では、近年開発が進んでいる指示文を使った並列プログラミング言語である OpenACC を使うことで、低開発コストでかつメ
ンテナンス性の高い、GPU 用大規模液状化解析手法を開発した。OpenACC を使うと、特定の構造の計算
であれば、CPU 用のコードにいくつかの指示文を挿入するだけで、GPU 上で並列計算を実装できる。このように、CPU 用のコードをベースにして、OpenACC を使って比較的容易に GPU 用コードを実装する
ことができるので、CUDA と比べて、開発コストが大幅に小さい。また、ほとんど CPU 用のコードと変わらないので、メンテナンス性が高い。つまり、OpenACC は実務への展開を達成するのに最適な選択肢の 1 つであるといえる。このように、比較的容易に GPU 用のコードを開発することが可能になってきているが、良い性能を達成するには、アルゴリズムについても考慮する必要がある。GPU のハードウェアの特性上、GPU を使うことによって高速化されやすい計算とあまり高速化されない計算があり、適切なアルゴリズムを使わないと、GPU の性能が発揮されないからである。一般に、本研究が対象とする非構造要素の有限要素法では、GPU が性能を発揮しづらい。それは、支配方程式の求解時に GPU が性能を発揮しづら
い計算の 1 つであるランダムメモリアクセスが多い計算を大量に行うからである。地震応答解析のような支配方程式が大規模疎行列方程式であるような場合の求解法には共役勾配法がよく用いられる。共役勾配法は反復解法の 1 つであるが、各反復において係数行列とあるベクトルの積を計算する必要がある。大規模疎行列とベクトルの積の計算では、ランダムメモリアクセスが大量に発生するため、GPU 計算の際に性能が低下する原因となってしまう。そこで、本研究では、OpenACC を使った大規模 3 次元地震応答解析
手法[5]を参考に、求解法に可変的前処理付き共役勾配法[9]を用いた(アルゴリズム 1: 7～17 行目)。可変的
前処理付き共役勾配法は、前処理で、前処理方程式を共役勾配法で解くという手法であり、通常の共役勾配法と比べ、反復回数が少なく[2]、ランダムメモリアクセスを抑えるため、GPU 計算に適した手法である。

65

図(2).18 構成則の計算アルゴリズム
構成則の計算についてもアルゴリズムを見直す必要がある。より現実に即した地震応答解析では、複数の材料が混在した複雑な地盤構造物モデルを解析する必要がある。そのため、構成則の計算時には、それぞれの要素について、どの構成則を使うのかを判定して、材料ごとに別々の計算を行うことになる。しかし、GPU は多数のコアが異なるデータに対して同じ演算を行う処理を行う際、優れた性能を発揮する一方、個々のデータに関して別々の演算を行う処理では、あまり性能を発揮できない。そのため、それぞれの要素について異なる計算を行う構成則の計算では、GPU の性能は発揮しづらい。そこで、CPU 上用に開発されたものではあるが、そのような条件分岐なしで、構成則の計算ができるような手法[10]を使用した。この手法の概要を図(2).18 に示す。この手法では、解析の最初にあらかじめ要素の順番を並び替えて置き、同じ構成則の計算をする要素が連続になるようにする。こうすることで、最初に材料 A の構成則計算を行い、次に材料 B の構成則を行い、最後に構成則 C の計算を行うというように、条件分岐なしで構成則の計算を行うことができる。この方法では、それぞれの材料ごとに、「異なるデータに対して同じ演算を行う処理」となり、GPU が性能を発揮できる。本研究の場合は、基盤層や構造物などの線形要素とマルチスプリングモデルに基づく非線形要素の 2 種類の構成則を使うため、プログラムの最初に、線形要素が配列の前半に、非線形要素が配列の後半になるように並び替えを行い(アルゴリズム 1: 2 行目)、構成則の計算の際は、まず、線形要素の計算(アルゴリズム 1: 20～25 行目)を行い、次に、非線形要素の計算(アルゴリズム 1: 26～40 行目)を行っている。このような要素の並び替えを行うことで、液状化解析での構成則計算も GPU による高速化が可能であり、前述の式(2-20)、(2-21)であらわされるような、300 本の 1 次元
ばねに関する計算を要素ごとに行うという膨大な計算コストを克服することができる。
GPU、CPU 上でのデータ管理も GPU 計算において重要である。GPU は CPU のメモリとは別にそれ自身のメモリを持っており、GPU 上で計算を行うには、GPU のメモリ上にデータが格納されている必要がある。一方で、MPI 並列などで並列計算を行う場合、他の MPI プロセスとの通信では基本的には CPU のメモリ上のデータをやり取りする。また、ファイルからの読み込みや、ファイルへの読み込みも CPU メモ
リを通して行う。そのため、CPU-GPU 間でデータの転送を行う必要がある。図(2).17(b)にハードウェア
間のデータ転送速度が示されている。CPU 上で計算を行うために、CPU メモリからデータを取得する場
合の転送速度は約 100 GB/s、GPU 上で計算を行うために、GPU メモリからデータを取得する場合の転送速度は約 1,000 GB/s であるのに対して、GPU-CPU 間でデータ転送を行う場合の転送速度は約 16 GB/s(倍精度実数 20 億個/s)、CPU 間で MPI 通信を行う場合の転送速度は約 10 GB/s とかなり遅い。GPU は 1 秒
間に最大約 8 兆回の倍精度実数演算を行うので、GPU-CPU 間のデータ転送を適切に行わないと、データ転送の時間が長くなり、計算速度が大幅に落ちる可能性がある。そこで、本研究では、ほとんどすべての計算を GPU 上で行うような実装を行った。この実装では、すべての変数をプログラムの最初に 1 度だけ
CPU から GPU にコピーし、それ以降は、基本的には、CPU-GPU 間のデータ転送は行う必要がない。これらのデータ管理は、OpenACC を使って実装できる。OpenACC では GPU 上での並列計算を指示すると、デフォルトでは、計算前に計算に必要なデータをすべて CPU から GPU に転送し、計算が終わった
66

ら、使ったデータすべての GPU から CPU に転送する。しかし、この中には、すでにこれまでの計算のた
めに、GPU メモリ上にコピーされていたデータや、GPU 上での計算の入力としては使われたが、値が変わっていないため、CPU メモリに戻す必要がないものなども含まれ、無駄なデータ転送が多い。本研究では、適切な OpenACC 指示文によって無駄なデータ転送を行わないような実装を行っている。
OpenACC を使ってデータ管理を行い、基本的には、すべての計算を GPU 上で行ってはいるが、並列計算時に他の GPU と通信が必要なデータやファイルへ出力する必要のあるデータなどは CPU-GPU 間で転送されている。この中でも、特に転送量が多いのが、支配方程式求解時の袖領域のデータである。有限要素法は陰解法であり、共役勾配法で、係数行列とあるベクトルの行列ベクトル積を並列計算する際(アルゴリズム 1: 15 行目、アルゴリズム 2: 7 行目)に、隣接する領域の袖領域の情報を通信する必要がある。そのため、その情報を送信側で GPU から CPU へ転送し、CPU 間で通信を行い、受信側で CPU から GPU へ転送する必要がある。この通信は、共役勾配法の各反復で必要であるため、各時間ステップで数百回程度の通信を行う必要がある。これらの CPU-GPU 間転送や CPU 間通信は、GPU 計算によって加速しないので、GPU 計算によって他の計算が加速されると、相対的に解析時間に占める割合が多くなり、高速化時のボトルネックになりうる。そこで、低精度変数 FP21[4]を通信に用いることで、転送量・通信量を減らし、この問題の克服を図った。本研究で採用している求解法である可変的前処理付き共役勾配法では、支配方程式の求解に至るまでに、アルゴリズム 1 の 9～17 行目 は 10 回程度反復されるが、前処理の計算(アルゴリズム 1: 11 行目)時の、Algorithm 2 の 5～13 行目での 1 時間ステップ当たりの合計反復回数は150～250 回程度である。つまり、袖領域のデータの通信が必要な行列ベクトル積計算のうち 94%～96%は、この前処理計算の中での行列ベクトル積計算(アルゴリズム 2: 7 行目)であるため、前処理での行列ベクトル積計算時の袖領域の通信を減らすことで、解析速度の向上が期待できる。この前処理計算は、反復回数の削減を図るためのものであり、解の精度には影響しないため、単精度(32bit)の変数による計算が行われている。
（求解のそれ以外の部分では、通常の科学計算で使われる倍精度変数(64 bit)が使われる）本研
究では、この前処理の行列ベクトル積計算(アルゴリズム 2: 7 行目)における袖領域の通信を 21 bit の低精度変数 FP21 に置き換えた。FP21 は、21 bit の実数変数で、通常使われる倍精度実数（64 bit）や単精度
実数(32 bit)よりも、精度が低く、データ量が少ない。前処理の計算は高い精度が必要ないので、FP21 を通信に使うことで、CPU-GPU 間のデータ転送と CPU 間の通信量を 2/3 に減らし GPU 計算におけるボト
ルネックになりうる通信時間の削減を図った。袖領域通信の際に送信側の GPU で単精度(32 bit)から低精度変数 FP21(21 bit)に変換し、FP21 で GPU から CPU にデータ転送し、CPU 間で通信を行い、受信側でCPU から GPU へデータ転送し、受信側の GPU で FP21 から単精度変数に戻して使用する。
③-(c). 数値解析
本研究で開発した GPU を使った大規模 3 次元液状化解析手法を使った数値実験を実施した。まず、開発手法の性能計測を実施し、CPU のみを使った液状化解析との解析時間の比較を行った。GPU を使うことで高速化が実現し、今まででは、現実的な時間で実行するには大規模な計算機環境が必要であった解析が、小規模計算機クラスタでも実施できるようになったことを確認した。適用例として、2 つの液状化対策について、それらの対策を模した地盤構造物モデルを使って液状化解析を行った。1 つ目の対策は、既設住宅の周りの地盤を格子状地盤改良によって強化するもの、2 つ目の対策は、河川の岸壁の背面地盤を液状化しない材料に置き換えるものである。本研究の開発手法によって、小規模な計算環境で、このような対策の性能照査を実施しうることが示された、解析対象の最大周波数は 2.5 Hz である。地表面の 3 次元形状が考慮されている。非線形層および液状化層は Iai の構成則[7]に従う。構造物と地盤、水はすべて四面体 2
次のソリッド要素でモデル化されており、水は、せん断剛性が 0 の弾性体としてモデル化されている。また、地盤、構造物、水の境界は結合している。このモデル化方法では、河川の応答がよく再現できるわけ
67

ではないが、河川が地盤や構造物に及ぼす圧力は十分再現できると期待される。今回の解析では、河川の応答が興味の対象ではなく、河川のからの圧力を受けている構造物・地盤がどのように応答するかが興味の対象であるので、このモデル化方法を使用した。
1) 計算時間の計測
我々は、CPU 上で本研究と同様の大規模 3 次元液状化解析を実施する検討を行ってきている[10]。そこで開発された手法で CPU のみを使って解析した場合と、本研究の開発手法で GPU も使って解析した場合の解析時間の比較を行った。図(2).19 に示した地盤モデルを使用した。表(2).3 に地盤物性を示す。モデルの自由度は 4,854,570、要素数は 1,179,817 である。地下水位は、地表面から 2m であり、入力波には、
1995 年兵庫県南部地震での観測波を使用し、時間ステップ幅 0.001 s で、100 time steps 分の液状化解析を行った。計算は、ABCI の 1 計算ノードを使用して行った。ABCI は GPU を搭載したスーパーコンピュ
ータで、各計算ノードに NVIDIA Tesla V100 (SXM2) GPU 4 台と、Intel Xeon Gold 6148 CPU 2 台が搭載されている。
表(2).3 性能計測用モデルの材料物性

図(2).19 性能計測用の地盤モデル
比較ケースは以下の 3 つである。CPU based: 先行研究[10]で開発された CPU 用の手法で、CPU のみを使用して解析、GPU based: 本研究の開発した GPU 用の手法（ただし、通信時に FP21 変数は使わない）で、GPU を使用して解析、GPU based+: 本研究の開発した GPU 用の手法（通信時に FP21 変数を使用）で、GPU を使用して解析。図(2).20 に各ケースの計算時間を示す。GPU based は CPU based と比較して、構成則の計算が約 14 倍、支配方程式求解が約 7.7 倍、全体として 9.9 倍の高速化を実現している。GPU based+ は GPU based と比較して、通信時間が 29%削減され、全体として 8%高速化している。GPU based+ は CPU based と比較すると 10.7 倍の高速化が実現した。

68

図(2).20 各ケースの計算時間後述する application 2 では 702m × 702m の領域を 最小 1m 間隔で離散化した 29,715,572 節点(89,146,716 DOF)の地盤構造物モデルに対して 30 秒間(時間ステップ幅 0.001s、30,000 time steps)の解析を行った。この解析に。既往研究[10]の手法を用いると、CPU-based のスーパーコンピュータ
Oakforest-PACS の 128 計算ノードを使って、14 時間 37 分かかっているが、本研究の手法を用いるとABCI の 13 計算ノードを使って 3 時間 33 分で実行できるようになっている。GPU を使って解析を高速化することにより、高詳細な 3 次元地震応答解析が小規模計算機クラスタ上でも実施できるようになることが示されている。
2) 適用例 1：格子状地盤改良による既設住宅の耐震性強化に関する解析
既存の住宅の周囲の地盤にモルタルなどを注入し格子状に地盤改良行うという液状化対策を模した解析を実施した。モデル概観を図(2).21 に示す。この地盤構造物モデルは住宅地を模したものである。モデルの自由度は 1,402,212、要素数は 331,944 である。地盤物性は表(2).4 に示す。建物の物性は、[11]を参考に、固有周期がおおむね建物の高さ × 0.1 s となるように設定した。建物の階数は建物の高さ/5 とした。
液状化のパラメータは性能計測時と同じである。地下水位は、地表面から 2m であり、入力波には、1995
年兵庫県南部地震での観測波[3]を使用し、時間ステップ幅 0.001 s で 30,000 time steps 分の液状化解析を行った。図(2).21(b)、(c)に示すように、一部の住宅に液状化対策を行い、対策の有無による結果の比較を通して、対策の有効性を確認した。GPU(tesla V100)4 枚を搭載した計算機クラスタで解析を行い、解析
時間は 1 時間 37 分だった。
表(2).4 適用例１の材料物性

69

(b) 地盤改良の位置(a) モデル概観(c)断面図図(2).21 適用例１の地盤構造物モデル
図(2).22 に Building A, B, C の変位の時刻歴応答を示す。Building A, B, C の第一固有周期はそれぞれ、0.90 s、0.35 s、0.29 s (x 方向)、0.64 s、0.78 s and 0.50 s (y 方向)である。Building A, B は対策をした住宅、Building C は未対策の住宅である。Building A, B において、対策によって、t=8 s ごろから、液状化によるものと思われる変位応答が抑制されている。一方で、対策をしていない建物 C では、変位はほとんど抑制されていない。対策によって、Building A, B での長周期成分(1.5 秒～)が小さくなり、スペクトルが全体的に短周期側に移動している。これは、対策によって液状化が抑制され、液状化によって引き起こされていた地盤の固有周期の長周期化が抑制されたからだと考えられる。x 方向に比べて y 方向のスペクトルで特に長周期成分が抑制されているのは、y 方向に対策を行った建物が連続しており、対策工全体
として、y 方向の抵抗が強くなっているからだと考えられる。また、未対策の建物 C においても、特に y方向の速度応答の長周期化がやや抑えられているが、これはおそらく隣接する建物 D の対策によって、y方向の応答が拘束されたからだと考えられる。図(2).23 は断面 EE'における、過剰間隙水圧比の最大値の分布を示している。対策によって、過剰間隙水圧の上昇が抑制されていることがわかる。この結果はこのよ
うな解析を 3 次元で行えるようになった意義が示されている。本手法を用いることで、格子状改良の位置や幅や材質を変えた際に、住宅の耐震性がどう変わるか、周りの住宅へのどのような影響を与えるのかな
ど、様々な検討が可能になりうる。

図(2).22 適用例２における Building A,B,C の変位応答

70図(2).23 適用例 2 における断面 EE’での過剰間隙水圧比の最大値の分布
3) 適用例 2：河川岸壁の背面地盤の置換
もう 1 つの液状化対策として岸壁の裏込め土の非液状化土への置換を模した解析を行った。
河口付近の実地形を模した地盤構造物モデルを使用した。モデルの概観を図(2).24 に示す。モデルの自由度は 89,146,716、要素数は 21,616,202 である。対策有の解析では、図(2).24(b)に示されるように、河川岸壁の裏込め土を非液状化土で置換した。地盤物性は 表(2).5 に示す。液状化のパラメータは性能計測時と同じである。入力波には、1995 年兵庫県南部地震での観測波を使用し、時間ステップ幅 0.001s、時間
ステップ数 30,000 で解析を行った。ABCI 13 ノードを使用して解析を行い、解析時間は 3 時間 33 分だった。
表(2).5 適用例 2 の材料物性(b) 拡大断面図(a)モデル概観

図(2).24 適用例２の地盤構造物モデル
岸壁の代表的な点として、岸壁の向きが異なる A, B, C での応答を図(2).25 に示す。対策によって、
t=13 s ごろからの液状化によると考えられる水平方向の変位応答が抑制されている。地表面における残留変位の分布を図(2).26 に示す。対策により、岸壁の残留変位が抑えられており、岸壁の耐震性が向上して71

いる。本手法を用いて、置換する土の物性や、置換する深さ・幅などを変えて解析することで、対策による耐震効果がどれくらいであるかを定量的に示すことができるようになりうる。

図(2).25 Point A, B, C における変位応答

図(2).26 地表面の残留変位の分布③-(d). まとめ開発された大規模数値解析コード群の実展開を考えると、大規模な 3 次元地震応答解析を小規模な計算機クラスタで実施できるようになることが必要である。本研究では、そのようなことが可能なハードウェ
アとして GPU に着目した。大規模な 3 次元地震応答解析の中でも複雑な問題として液状化解析を例に、OpenACC を使って、GPU 用の解析手法を開発し、比較的容易に 10 倍もの高速化を実現した。これにより、先行研究では、CPU ベースのスーパーコンピュータ上で 128 計算ノードを使って 14 時間 37 分かかっていた解析が、GPU を搭載したスーパーコンピュータの 13 計算ノードを使って 3 時間 33 分で実行できるようになり、大規模な 3 次元応答解析を小規模な計算クラスタで実施できるようになった。また、大規72

模 3 次元液状化解析が小規模な計算クラスタで実施できるようになったことにより、可能になるであろう
検討の例として、2 つの液状化対策を模した解析をおこなった。小規模な計算機で大規模な 3 次元液状化解析が可能になったことで、液状化対策の様々な性能照査が容易になり、対策の選定や、対策法の開発、改良などに役立つと考えられる。
参考文献[1] Tsuyoshi Ichimura, Kohei Fujita, Seizo Tanaka, Muneo Hori, Maddegedara Lalith, Yoshihisa Shizawa
and Hiroshi Kobayashi, Physics-based urban earthquake simulation enhanced by 10.7 BlnDOF x 30 K
time-step unstructured FE non-linear seismic wave simulation, SC14: International Conference forHigh Performance Computing, Networking, Storage and Analysis, pp 15-26, 2014. doi:
10.1109/SC.2014.7[2] Tsuyoshi Ichimura, Kohei Fujita, Pher Errol Balde Quinay, Lalith Maddegedara, Muneo Hori, Seizo
Tanaka, Yoshihisa Shizawa, Hiroshi Kobayashi and Kazuo Minami, Implicit Nonlinear Wave
Simulation with 1.08T DOF and 0.270T Unstructured Finite Elements to Enhance Comprehensive
Earthquake Simulation, SC15: International Conference for High Performance Computing,
Networking, Storage and Analysis.
[3] Tsuyoshi Ichimura, Kohei Fujita, Masashi Horikoshi, Larry Meadows, Kengo Nakajima, Takuma
Yamaguchi, Kentaro Koyama, Hikaru Inoue, Akira Naruse, Keisuke Katsushima, Muneo Hori,
Maddegedara Lalith, A Fast Scalable Implicit Solver with Concentrated Computation for
NonlinearTime-evolution Problems on Low-order Unstructured Finite Elements, 32nd IEEE
International Parallel and Distributed Processing Symposium, 2018.
[4] Tsuyoshi Ichimura, Kohei Fujita, Takuma Yamaguchi, Akira Naruse, Jack C. Wells, Thomas C.
Schulthess, Tjerk P. Straatsma, Christopher J. Zimmer, Maxime Martinasso, Kengo Nakajima, Muneo
Hori, Lalith Maddegedara, A Fast Scalable Implicit Solver for Nonlinear Time-Evolution Earthquake
City Problem on Low-Ordered Unstructured Finite Elements with Artificial Intelligence andTransprecision Computing, SC '18 Proceedings of the International Conference for High Performance
Computing, Networking, Storage, and Analysis, 49:1-49:11, 2018.
[5] Yamaguchi T, Fujita K, Ichimura T, Naruse A, Lalith M, Hori M. GPU implementation of a
sophisticated implicit low-order finite element solver with FP21-32-64 computation using OpenACC.
In: Proceedings of Sixth Workshop on Accelerator Programming Using Directives. 2019.
[6] Iai S, Matsunaga Y, Kameoka T. Strain space plasticity model for cyclic mobility. Soils andFoundations 1992; 32(2): 1-15.
[7] Iai S. Three dimensional formulation and objectivity of a strain space multiple mechanism model forsand. Soils and Foundations 1993; 33(1): 192-199.
[8] Mu D, Chen P, Wang L. Accelerating the discontinuous Galerkin method for seismic wave propagation
simulations using multiple GPUs with CUDA and MPI. Earthquake Science 2013; 26: 377–393.
[9] Golub GH, Ye Q. Inexact preconditioned conjugate gradient method with inner-outer iteration. SIAM
Journal on Scientific Computing 1999; 21(4): 1305-1320.
[10] Kusakabe R, Fujita K, Ichimura T, Hori M, Wijerathne L. A Fast 3D Finite-Element Solver forLarge-Scale Seismic Soil Liquefaction Analysis. In: Computational Science – ICCS 2019. Springer
International Publishing; 2019; Cham: 349–362.
73

[11] Isbiliroglu Y, Taborda R, Bielak J. Coupled Soil-Structure Interaction Effects of Building Clusters
during Earthquakes. Earthquake Spectra 2015; 31(1): 463-500.
（３）開発した計算アプリケーション群を実務で活用するための環境整備（１）において整備している国の被害想定算出に使われるものと同等の長周期地震動計算アプリケーシ
ョン群の実務への活用を進めるため、建築・土木系の各企業が独自の問題設定で実施するためのサポート
体制を、海洋研究開発機構が中心となって整備している。本年度は、本課題の連携企業の皆様が計算アプ
リケーション群を使用するための環境構築を進めた。計算アプリケーション群の基幹数値解析コンポーネントである E-wave FEM の計算環境を、連携企業の皆様が利用可能な形態で、構築した。更に、当該計算環境を有効に利用して頂くため、連携企業の皆様を対象に、5 回からなるウェビナー形式（Web 上のオン
ラインセミナー）のチュートリアルを実施した。その上で、連携企業の皆様から E-wave FEM 計算環境に関するフィードバックを頂いた。本稿①にて、これらの取り組みについて述べる。また、連携企業の皆様には、企業でニーズのある地震動計算への活用を進めていただいている。その中ですでに得られた成果について②で述べる。
① 実施した環境整備の概要①-(a). E-wave FEM の計算環境の構築
計算アプリケーション群の基幹数値解析コンポーネントである E-wave FEM の計算環境として、大規模共有メモリ計算機及び高並列スーパーコンピュータからなる計算システムを構築した（図(3).1、図(3).2）。
この計算システムは、海洋研究開発機構が所有し理化学研究所計算科学研究センターに設置された大規模共有メモリ計算機”Jasper”と、東京大学情報基盤センターと筑波大学計算科学研究センターが共同運営する、最先端共同 HPC 基盤施設（JCAHPC：Joint Center for Advanced High Performance Computing）の共同利用（高並列）スーパーコンピュータ”Oakforest-PACS”(https://www.cc.utokyo.ac.jp/supercomputer/ofp/system.php)とを、主構成要素とする。エンドユーザは、インターネットを経由した SSH 接続（共通鍵暗号方式）により、この計算システムを利用する。
①-(b). マニュアルの整備
計算環境の構築に合わせ、E-wave FEM のマニュアルを作成した。作成したのは、「基本マニュアル」と「チュートリアルマニュアル」の二種類である。基本マニュアルは、各プロセスと入力用計算設定ファイ
ル・出力ファイルとの対応、及び個別ファイルの内容を詳細に説明している。但し、その説明は総花的な羅列を避け、利用手順に沿った実用的な記述としている。今後予想される計算機の利用環境を含む計算環
境の変更、計算プログラムの改良等に短時間で対応するため、基本マニュアルは頻繁にアップデートし、連携企業の皆様に提供を続ける。チュートリアルマニュアルは、①-(c)で説明する全 5 回のチュートリアルにて実施された解説のベースとなる資料である。ウェビナー録画とともに、協賛企業参加者の皆様がアクセス可能なサイトに公開されている。
①-(c). チュートリアルの実施
チュートリアルの概要を表(3).1 に示す。2020 年 10 月から 12 月までの期間に、本課題の連携企業の皆様を対象に、本課題にて開発された地震動計算用大規模有限要素計算プログラム「E-wave FEM」の使用方法についてのチュートリアルを実施した。チュートリアルは全 5 回により構成され、有限要素計算用メッシュの生成、並列計算用の領域分割、スーパーコンピュータ上での高並列地震動計算、計算結果の可視化までをウェビナー形式で説明を実施した。一連のチュートリアルでは、説明と同時に、参加者の皆様に74

はハンズオン（遠隔実習）により、実際に大規模共有メモリ計算機”Jasper”、及び高並列スーパーコンピ
ュータ”Oakforest-PACS”に SSH 接続し、プリプロセス、E-wave FEM による地震動計算、ポストプロセ
スの一連の流れ(図(3).2)を経験して頂いた。
①-(d). 計算環境に対するフィードバック
チュートリアルの参加者から頂いたフィードバックの概要を表(3).2 に示す。一連のチュートリアルを通じて連携企業の皆様に E-wave FEM をご使用いただき、その後複数回のオンライン会議、メール等によりE-wave FEM 計算環境に関するフィードバックを頂いた。フィードバックは、チュートリアルによる Ewave FEM 計算環境に関する設定、手順等の直接的な使用感に関するものと、オンライン会議等で議論された E-wave FEM 計算環境の今後の利用に関するものがあった。ここでは、それら二つを順に紹介する。
E-wave FEM 計算環境に関する設定、手順等に関するフィードバックは下記のとおりである：


E-wave FEM による地震動計算は、3 プロセス（プリプロセス、地震動の有限要素計算、及びポストプロセス）から構成されている。各プロセスは機能毎にモジュールとして分割され、そのモジュール
をプロセス毎に連続処理する構成となっている。連携企業の皆様からは、機能毎のモジュール化とそ。
の連続処理によるエンドユーザの負荷軽減とに関し、評価するコメントを頂いた（＊1）



一方、改善すべき数々の問題点も指摘して頂いた（＊2）。それらの問題点は、大きく分けて 5 つにまとめられる。第 1 にミスをなくす機構について、第 2 に計算設定（多数回計算、減衰）について、第3 に震源設定（特に面震源と震源過程）について、第 4 に地盤モデル設定（高周波地震動計算と地盤モデルとの関係、メッシュの自動生成）について、及び第 5 に計算結果の評価について、である。



第 1 のミスをなくす機構に関しては、各モジュールで処理がどの様に進行しているかを示すログ、及び画面上への警告表示が不十分との指摘がされた。各モジュールはその処理毎にログを記録しているが、どのログが何を記録しているかの説明（チュートリアルでの）が不足しているとの指摘を受けた。



第 2 の計算設定に関しては、多数回計算を前提としたインターフェイスの必要性、減衰設定の自由度の低さについて指摘を受けた。特に、減衰に関し、現在実務で広く利用されている有限差分法との減衰設定の違いを指摘するコメントを頂いた。



第 3 の震源設定に関しては、面震源要素の変更と地盤モデル再構築との関係と、震源時間関数設定に関し指摘を受けた。また、内陸の活断層（例：高角、横ずれ成分が優勢）を想定した断層モデルの設定を期待するコメントを頂いた。



第 4 の地盤モデルに関しては、高周波地震動計算と地盤モデル精度との関係、全自動メッシュ生成ツールの有効性について、地盤モデルへの不均質性の導入について、のコメントを頂いた。



第 5 の計算結果の評価に関し、「検証と妥当性確認（Verification and Validation）」リポートを期待するコメントを頂いた。



これらの改善点に関するコメントは、E-wave FEM 計算環境を更に改善する上で重要な指摘であり、指摘された改善点の幾つかは既に実装を開始している。
今後の利用に関して（＊3）は、複数回のオンライン会議により、協賛企業の皆様から今後への改善点、ご提案を含めたフィードバックを頂いている。一部前述のフィードバックと重なるが、改めてそれらを紹介する：


震源から地震基盤までの地震動を本課題で開発、運用される E-wave FEM 計算環境により計算し、その地震動を利用して工学的基盤から上を協賛企業各社が独自に解析する、という使い分けがいいので
はないか。

75



この協賛企業各社の共通基盤となる E-wave FEM による地震動計算は、本課題の最新の計算手法、地盤データ等によってその信頼性が担保されることが望ましい。



その信頼性を維持するためにも、逆問題解析による地盤モデルの更新、計算プログラムの改良等を通して、共通基盤部分のアップデートをしてほしい。



海水を考慮した地盤モデルで計算できないか。また、地盤モデルに海洋研究開発機構が収集解析した
海底地形データを使えないか。



今後の研究的な側面として、この計算システムを三次元不均質構造でのグリーン関数計算に利用した
い、あるいは震源インバージョン等の逆問題解析の構成要素として利用したい。



スーパーコンピュータ富岳との連携を考慮して、本計算システムのインターフェイスを富岳のインターフェイスと揃えてはどうか。
これら改善点のご指摘、今後の利用に関するご意見、ご提案を踏まえて、連携企業の皆様との議論を深めながら、この計算システムを社会実装に足る機能とインターフェイスを備えたアプリケーションとして整備を進めていく。

図(3).1 E-wave FEM 計算環境

図(3).2 E-wave FEM 計算の流れ

76表(3).1 チュートリアル実施内容
日付・形式

内容

第1回

（E-wave FEM と計算環境の紹介・水平成層構造モデル作成）

2020 年 10 月 8 日１）E-wave FEM の概要

ウェビナー形式

２）事前準備・ 計算プロセスの紹介
・ 大規模共有メモリ計算機(Jasper)、高並列スーパーコンピュータ(Oakforest-PACS)紹介
・ 各計算機への接続設定
３）地震動計算のプリプロセス（Jasper 上）・ 水平成層構造モデルの紹介
・ 水平成層構造モデルの形状・物性の設定
・ 水平成層構造モデルのメッシュ作成
・ 高並列計算用地盤モデルの分割
・ 分割された地盤モデルを Oakforest-PACS へ転送

第2回（計算環境整備の確認・複雑構造モデル(Moon Basin)の作成）

2020 年 10 月 20 日

１）事前準備（第 1 回のレビュー）

ウェビナー形式

２）地震動計算のプリプロセス（Jasper 上）・ 複雑構造モデル(Moon Basin)の紹介
・ 複雑構造モデル(Moon Basin)の形状・物性の設定
・ 複雑地盤構造モデルのメッシュ作成
・ 高並列計算用地盤モデルの分割
・ 分割された地盤モデルを Oakforest-PACS へ転送

第3回（Oakforest-PACS 上の地震動計算（点震源）とポストプロセス）

2020 年 11 月 10 日

１） プリプロセス（水平成層・複雑構造モデル作成）のレビュー

ウェビナー形式

２） 点震源を含む問題設定の説明
３） 水平成層構造モデルと点震源による地震動計算
・ 水平成層構造モデル上の計算環境、点震源における震源過程、及び観測点の設定
・ 地震動計算の実行
・ ポストプロセスの実行・観測点における波形図の生成
３）複雑構造モデルと点震源による地震動計算
・ 複雑構造モデル上の計算環境、点震源における震源過程、及び観測点の設定
・ 地震動計算の実行
・ ポストプロセスの実行・観測点における波形図の生成
４）結果の取りまとめ

第4回（Oakforest-PACS 上の地震動計算（面震源）とポストプロセス）

2020 年 11 月 18 日

１）第 1 回－第 3 回までのレビュー

ウェビナー形式・ 計算環境・接続設定の確認
・ プリプロセス・本解析（地震動計算）・ポストプロセスの流れ確認
77

・ 水平成層構造モデル・複雑構造モデル各々によるメッシュ生成、分割、転送の確認
２）複雑構造モデルと面震源による地震動計算
・ 面震源を含む問題設定の説明
・ Jasper 上の計算環境、面震源に対応した複雑構造モデルの設定
・ 複雑構造モデルの生成・分割と Oakforest-PACS への転送
・ 複雑構造モデル上の計算環境、面震源における震源過程、及び観測点の設定
・ 地震動計算の実行
・ ポストプロセスの実行・観測点における波形図の生成
３）結果の取りまとめ
第5回（地震動計算流れ確認、面震源設定、可視化ツール導入、質問への回答）

2020 年 12 月 15 日

１） プリプロセス、地震動計算、ポストプロセスの確認

ウェビナー形式

２） 面震源設定の補足説明
・ 面震源の詳細説明
・ 面震源用データ作成プログラムの説明
３） 計算波形の可視化ツール導入
４） 質問への回答

表(3).2 フィードバック（コメント概要）
評価するコメント

１）計算システムの構成（＊１）

・ 「大規模計算環境とそれを使えるプログラムの組み合わせ」というところに強みがある。
２）計算システムの手続き、使用感・ 計算手順はわかりやすかった。一括実行シェルも使いやすい。
・ メッシュ生成が全自動で使いやすい。
・ 使用感としては差分法と比べてそんなに変わらない。
・ プログラムは分かりやすかった。うまく構成されていた。人の手をあまり入れずにできた。
・ 地盤モデルの層境界や震源時間関数をサブルーチンで与えるのは自由度が高くて良い。
・ 可視化部分もやりやすい。従来の差分計算系のプログラムと遜色がない。

改善点を指摘する1.

コメント（＊２）・ ログ・警告が少ない。ログファイルの見方のガイダンスが足りない。

ミスをなくすための機構

・ 一目で設定の概要がわかるファイルが欲しい。
・ 所要時間のわかるログを表示して欲しい。
・ メッシュサイズの自動調整の際の間違えやすさ（警告だけでよい。少なくともメッシュサイズを切替え可能に。
）
2.

計算設定

・ 多数回計算78長周期に対しても、大量に計算を回してわかることがあるはず。
回数はターゲット、パラメータ数次第。
多数回計算のしやすいインターフェイスの必要性がある。
・ 計算の待ち時間が長い。
・ 減衰のより柔軟な設定（例えば質量比例減衰，S 波と P 波で別々の減衰など）ができると汎用性が高まる。
・ 周波数範囲（FMAX, FMIN）, Q 値を入力したときに、どのような振動数特性を持つ減衰になるのかというレイリー減衰の説明が記載されていると親切だと思う。
・ 地盤モデルを生成した後でも、解析時にレイリー減衰を設定できるようにし
てほしい。
3.

震源設定

・ 破壊伝搬速度等、研究のための細かい設定をする部分に、あまり自由度がな
い。
・ 断層面の変更にはメッシュの切り直しが必要になる。
・ できれば内陸の活断層（高角の断層面の設定）にも使えるとよい。
・ 震源時間関数がサブルーチンとして予め準備されていない形状であっても自力で入力できた。今後は正式にサポートしてほしい。
4.

地盤モデル設定・ 高周波までの計算を考慮すると、深部地盤データの精度が足りない。
設計用地震動（特に高周波）はどの様に計算するのか。
震源はどんな高周波を励起するのか。
深部構造はどんな高周波を励起するのか。
・ 全自動メッシュ生成ツールは、地盤モデルの境界が入り組んだ時にどの様になるかは使ってみないとわからない。
・ 複雑な形状を扱うとよりツールのことが分かってくると思う。
・ 地盤モデルに層内不均質・ランダム不均質を導入できるとよい。
・ 数式で表せない層境界の外部入力ができるとよい。
5.

計算結果の評価

・ V&V のレポートがあると重要構造物の解析で使いやすくなる。
今後の利用へのコ
メント（＊３）・ 震源から理学的基盤面までを想定地震動として E-wave FEM によって計算し、基盤面より上側を各社が独自に解析する、という切り分けがよいのではないか。
・ 三次元不均質構造でのグリーン関数計算、震源インバージョンに使いたい。
・ 海域（海水層を入れた）での計算で、JAMSTEC の探査した海底地形は使用できるか。
・ 津波の計算はできるのか。
・ 今後富岳を使う必要性が高まった場合、富岳を利用した計算と、OakforestPACS を利用した計算とのインターフェイスが予めそろっているとよい。
79

② 企業においてニーズのある地震動計算への適用―モンテカルロシミュレーションを用いた地震動評価―
チュートリアルにより E-wave FEM の使い方を習得された連携企業の皆様に、企業でニーズのある地震
動計算への活用を進めていただいている。その中ですでに得られた成果について述べる。
近年、確率論的地震ハザード評価の精緻化を目的として、断層モデルを用いた地震動評価の導入が検討されている。[1]は、熊本地震前震（2016 年 4 月 14 日、Mj6.5）について、モンテカルロシミュレーション（MCS）により震源特性の不確定性を考慮した地震動を作成し、観測記録との比較や微視的震源特性が短周期の地震動に及ぼす影響について分析している。長周期地震動についても、今後発生が懸念される南海トラ
フや相模トラフの巨大地震等に対して、震源特性のばらつきを考慮した地震動評価が重要となる。3 次元不均質地盤構造における地震動を精度よく、大規模計算機において高速に実行できる E-wave FEM は、ここに存在する現実的な地盤構造を導入した長周期地震動計算を多数実施するニーズに合致している。ここでは、[1]による熊本地震前震の震源特性を用いて、長周期地震動の地震波群を E-wave FEM により作成し、震源特性の不確定性が長周期地震動に及ぼす影響について考察した。
地下構造モデルは、地震基盤までは J-SHIS の深部地盤構造モデルを、地震基盤以深は地震調査研究推進
本部の全国 1 次地下構造モデルの地殻構造を用いた。図(3).3 にモデル化範囲、断層位置、評価地点位置を、
表(3).1 に使用した地盤物性を示す。解析周期を 2.5 秒以上として 1 波長当たり最低 5 要素の計算条件でメッシュ生成したところ、最小要素サイズは 250m となり、約 2 億 1160 万節点、約 1 億 5730 万要素の計算モデルが得られた。作成した計算モデルを、並列計算ができるよう 1024 個に領域分割した。時間ステッ
プ幅を 0.2 秒として、410 ステップ、つまり 82 秒分の計算を行うこととした。
評価地点は、K-NET、KiK-net の観測点のうち観測記録が得られている地点とした。[1]では、巨視的震源特性は[2]を用い、微視的震源特性に対しては MCS で不確定性を与え、統計的グリーン関数法を用いた短周期成分を対象とした地震動評価を行っている。本検討でも、これと同じ震源特性を用い、長周期成分を対象とした地震波群を作成した。表(3).2 に、巨視的震源特性を、表(3).3 に、[1]による微視的震源特性を示す。
図(3).4 に、断層モデルと震源特性の対応図を示す。その他の震源特性は[3]に従って設定した。本検討では、
背景領域は考慮せず、アスペリティのみの評価とした。また震源時間関数は二等辺三角形とした。以上の震源特性の不確定性を考慮した MCS を実施し、各評価地点について 1000 波の合成波形群を作成した。
このような計算・震源モデルを用いて長周期地震動計算を実施し、Oakforest-PACS の 128 ノードを計
168 時間用いてすべての計算を完了した。図(3).5 に、KMM001、KMM018、MYZH04 の結果を示す。図(3).3 に示す通り、KMM001 と KMM018 は断層走向方向の延長線上に位置し、MYZH04 は断層走向と直交する方向に位置する。震央距離は約 49～55km である。図中の灰色線は合成波形群の擬似速度応答スペクト
ル(h=5%、以降 pSv)、黒線がその中央値、赤線が観測記録の pSv である。図(3).6 に合成波形群の各周期における自然対数標準偏差(S.D.)を示す。本検討の S.D.は、[1]の短周期と比較して大きい傾向がみられる。図(3).7 に中央値と S.D.から算出される観測記録の偏差を示す。観測記録の偏差は合成波形群の±3σ の範囲に収まっており、MCS で考慮した不確定性の条件はやや遠方の観測記録の振幅を説明可能であると考えられる。ただし、近傍の観測点ではより複雑な震源特性の影響を受けるため、合成波形群と観測記録が乖離する
観測点も多く、記録との整合性に関しては、今後の検討が必要である。
線形回帰分析を実施し、各震源特性が合成波形へ及ぼす影響について確認した。各変数を対等に評価する
ため、平均値 0、標準偏差 1 の正規分布に正規化した上で実施した。図(3).8 に各震源特性の回帰係数を示す。震源特性のうち、破壊開始点 startX、アスペリティ位置 aspX1、アスペリティすべり量比 Cda が高い感度を示すことが確認できる。また startX、aspX1 は観測点により差があるが、CVr、Cda、Cτa は断層との位置関係の影響をあまり受けていない。
破壊開始点、アスペリティと観測点の相対的位置関係について、
KMM001は走向反対方向に位置することから、startX が正の値で増加する場合には、破壊が KMM001 に向かって進展する。また aspX1 が正の値で増加する場合には、アスペリティが KMM001 から遠ざかる。そのため、
80startX は正の相関、aspX1 は負の相関になっていると考えられる。KMM018 は走向方向に位置し、KMM001と逆の結果となり、MYZH04 は走向直交方向に位置し、startX の相関は小さい。断層の破壊進展方向が大きく影響していることが確認できる。アスペリティすべり量比 Cda について、アスペリティのすべり量はア
スペリティの地震モーメントに比例することから、地震規模も長周期地震動に対して及ぼす影響が大きいと考えられる。本検討を[1]の短周期の検討と比較すると、理論あるいは数値解析的な地震動評価は、相対的位
置関係や震源パラメータの変動の影響を受けやすいと考えられ、評価結果のばらつきが大きくなる傾向がある。特に、破壊進展方向の影響が大きい傾向がある。その他の特性については、調和的な結果となった。
熊本地震前震の観測記録と MCS による長周期地震動の合成波形群を比較し、震源特性の不確定性に起因するばらつきによりやや遠方の観測記録の振幅レベルを説明できることを確認した。また、MCS の結果に対し線形回帰分析を行い、短周期に比べ、相対的位置関係や地震規模といった震源特性が長周期地震動に与える影響が大きいことを確認した。
このように、地震動計算をモンテカルロシミュレーションに適用することにより、長周期地震動に影響を与える震源特性に関する知見を得ることができ、大規模計算機において地震動計算を高速に行える E-wave
FEM の有用性が示されたといえる。

モデル化範囲KMM001
断層基準点
震源（観測）

巨視的断層面

MYZH04
KMM018

図(3).3 モデル化範囲、断層位置、評価地点の位置

81

表(3).1 地盤物性
No.
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

Vp[km/s]
1.80
1.80
1.80
1.80
1.80
2.00
2.00
2.10
2.10
2.20
2.30
2.40
2.40
2.50
2.50
2.60
2.70
3.00
3.20
3.40
3.50
3.60
3.70
3.80
4.00
4.00
5.00
4.60
5.00
5.50
5.50
5.70
6.00
6.40
7.50
5.00
6.80
8.00

Vs[km/s]
0.50
0.50
0.50
0.50
0.55
0.60
0.65
0.70
0.75
0.80
0.85
0.90
0.95
1.00
1.10
1.20
1.30
1.40
1.50
1.60
1.70
1.80
1.90
2.00
2.10
2.10
2.70
2.90
2.70
3.10
3.20
3.30
3.40
3.80
4.50
2.90
4.00
4.70

ρ[g/cm3]
1.90
1.90
1.90
1.90
1.90
1.90
1.95
2.00
2.00
2.00
2.05
2.05
2.10
2.10
2.15
2.15
2.20
2.25
2.25
2.30
2.30
2.35
2.35
2.40
2.40
2.40
2.50
2.55
2.50
2.60
2.65
2.70
2.75
2.80
3.20
2.40
2.90
3.20

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
300
300
300
400
500
200
300
500

表(3).2 巨視的震源特性
震源位置32.7417、130.8087

震源深さ11.39 km

地震モーメント

2.04×1018 Nm

走向、傾斜角212°/ 89°

すべり角

−164°

断層長さ

14 km

断層幅

13 km

82

参照

J-SHIS全国 1 次地下構造モデル

表(3).3 微視的震源特性
震源特性

記号

分布形状

λ (μ)※1

ζ (σ) ※1

破壊伝播速度のせん断波速度比

CVr

対数正規

0.72

0.08

ライズタイム係数

Cτa

対数正規

0.50

0.20

Csa

正規

(0.22)(0.04)

アスペリティ面積比

※4

Nasp

アスペリティ個数
アスペリティの面積比※2

アスペリティすべり量比

Csa12

正規

(0.50)(0.30)

Cda

正規(2.00)(0.68)

startX一様

アスペリティ下端に配置

一様

断層面内に配置※3

破壊開始点

断層内に 1 個または 2 個

aspX1、aspY1

アスペリティ位置

aspX1、apsY2※1 λ：中央値、μ：平均値、ζ：自然対数標準偏差、σ：標準偏差※2 いずれかのアスペリティ面積が 10％を下回る場合にはアスペリティ個数は 1 個とした
※3 アスペリティのすべり量が背景領域のすべり量を下回らないように下限値を設定※4 背景領域の地震モーメントが負にならないよう Csa と Cda の条件を設定
北
断層基準点

走向方向

走向θ
aspY2
アスペリティ2

傾斜角δ

aspX1

aspX2
背景領域

アスペリティ1

aspX1
startX

巨視的断層面

破壊開始点

図(3).4 断層モデルと震源特性の対応図

83

EW

NS

UD

（a）KMM001
EW

NS

UD

（b）KMM018
EW

NS

UD

(c) MYZH04
図(3).5 合成波形の応答スペクトル

EW

NS

KMM001

UD
KMM018

MYZH04

図(3).6 合成波形群の自然対数標準偏差

EW

NS

KMM001

UD
KMM018

図(3).7 観測記録の偏差84

MYZH04

startX

aspX1
aspX2

aspY1
aspY2

Csa
Csa12

Cτa
EW

Cda
NS

CVr
UD

KMM001

KMM018

MYZH04

KMM001

KMM018

MYZH04

図(3).8 震源特性の地震動応答に対する感度（回帰係数）
参考文献[1] 青木雅嗣、五十嵐さやか、内山泰生、坂本成弘：モンテカルロシミュレーションにより不確定性を考慮した断層モデル地震波作成手法の妥当性検証 －2016 年 4 月 14 日熊本地震前震(Mj6.5)の例－、日本建築学会大会学術講演梗概集（東北）、構造 1、pp.85-86、2018[2] Asano and Iwata：Source rupture processes of the foreshock and mainshock in the 2016 Kumamoto
earthquake sequence estimated from the kinematic waveform inversion of strong motion data、Earth,Planet Space 68、Article number: 147、2016[3] 地震調査研究推進本部：震源断層を特定した地震の強震動予測手法 (「レシピ」)、2020（４）プロジェクトの総合的推進
プロジェクトの総合的推進のために、プロジェクトに参加する研究者全体の連携を密とすることは必須であるが、運営に過度の労力を割くことは決して得策ではない。そこで、全体の効率的かつ円滑な運営のため、まず始めに、実施機関に加え連携機関である内閣府(防災担当)や建設関係企業、大学・国立研究開発
法人を含めた参画機関全体でのキックオフ全体会議(2020 年 6 月 11 日 Web 形式)を開催し、課題の目標、計画などの情報を代表機関より共有した上で、メンバー紹介やプロジェクトへの期待・要望に関する意見交換を行った。その後、実施機関での進捗管理等のための会合を定期的に開催した（隔月開催の全体会
議）。さらに、外部有識者を招聘した諮問委員会を開催した（半期に一度開催）。
プロジェクト全体の協力・連携を密とし円滑に運営していくことで、より効果的に成果を創出できるよう下記の諮問委員会を組織化した。諮問委員会からは、本プロジェクトの進捗状況の把握・評価・改善提言・指導等を受けた。

85

委員長
平田 直

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

委員
渦岡 良介

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

小野 謙二

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

山中 浩明

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

矢来 博司

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

委員長・委員の専門は、観測地震学・地盤工学・地震工学・数値流体力学・計算科学・物理探査工学・
強震動地震学・固体地球科学である。また、本課題で構築される統合的予測システムを社会実装する際、最大のユーザとなることが想定される国レベルの地震防災・減災の有識者でもある。
第 1 回の諮問委員会は 2020 年 10 月 2 日に Web 形式で開催された。この委員会では、プロジェクト概
要および研究計画の説明と約半期での進捗状況を説明した。本年度が本プロジェクト実施の初年度である
ため、代表機関と協力機関其々の実施機関の事業進捗状況の報告に対する把握、提言・指導が行われた。
実施計画に関しては様々な意見が出されたが、好意的な意見が多く、概ね、研究計画の妥当性・合理性・重要性が了解され、本プロジェクトへの期待が寄せられた。第 2 回の諮問委員会は 2021 年 3 月 8 日にWeb 形式で開催された。この委員会では、本年度の研究成果を報告し、評価を受けた。特に本プロジェク
トにて推進する、SC20 での発表内容を中心とした計算科学としての新規性については高く評価された。また、社会実装を目指す観点から、国民目線での成果の出し方について助言を受けた。具体的には、「今まで計算能力の関係で簡便化した手法・モデルでやっていたものを、現実のモデルに近づけるというところが非常に重要。そこを（高性能計算技術を使うことで）胸を張って堂々と是非やっていただきたいと思う」「その上で、それだけ精緻なハザードの予測を社会がどう利用するかところで、関係プロジェクトとの連携の議論が重要となる」という趣旨のコメントを受けた。結論として、本年度の研究の進捗状況に関して、特段の問題がないことが了解された。
プロジェクトで得られた成果については、広報等を通じて積極的に公表し、今後の展開に資することを進めた。具体的には本プロジェクトの和文・英文ホームページ（図(4).1）を介した公表である。イベン
ト・研究成果について最新情報を掲載するページを設け、2021 年 4 月までに合計 20 点の記事を記載した。また、本年度の成果である「Green's function library for subduction zones」についても、Web 上での公開をおこなっている。さらに、令和２年度成果報告会(2021 年 3 月 8 日 Web 形式)を開催した。成果
報告会では本年度の成果に係る報告の後、連携企業も交え闊達な議論がおこなわれ、地震に関する災害被害予測アプリケーションの実用化に向けての提案もあった。
計算結果の可視化は広報に重要である点を重視した。前身である重点課題より引続き、理化学研究所計算科学研究センターに設置したサーバ一式を関連するポスト処理に利用することを継続した。実都市を対象とするため、不特定多数が利用するホームページには公表できないものの、一般を対象としたシンポジ
ウムを含む会議等で可視化結果の動画を公表した。

86

図(4).1 和文・英文ホームページイメージ
計算機環境については、理化学研究所計算科学研究センターおよび、海洋研究開発機構に整備した計算機環境により、重点課題で開発したアプリケーション群や各関係機関で得られた貴重なデータならびに、本年度の業務で得られた成果等の保管・管理がなされている。また、研究を推進するための計算機環境の維持・整備として、これら計算機の定常的な監視や定期保守を実施した。

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

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

開催日
令和 2 年 6 月 11 日
令和 2 年 7 月 17 日
令和 2 年 9 月 4 日
令和 2 年 11 月 6 日
令和 3 年 1 月 15 日
令和 3 年 3 月 3 日

委員会

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

令和 2 年 10 月 2 日
令和 3 年 3 月 8 日

成果報告会

令和 2 年度成果報告会（オンラ

令和 3 年 3 月 8 日

イン）
ワークショップ

令和 2 年 10 月 8 日

チュートリアル令和 2 年 10 月 20 日
令和 2 年 11 月 10 日
令和 2 年 11 月 18 日
令和 2 年 12 月 15 日

87

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

担当機関

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

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

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

堀

アプリケーション活用

海域地震火山部門

研究員

担当責任者

縣

高峰
亮一郎

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

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

教授

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

東京大学地震研究所

市村 強ーションの開発

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

准教授
藤田 航平

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

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

センター長

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

堀備

海域地震火山部門

研究員
縣

高峰
亮一郎

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

センター長

推進

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

堀

海域地震火山部門

88

高峰

別添１

学会等発表実績

1．学会等における口頭・ポスター発表
発表した成果

発表した場所
発表者氏名（発表題目、口頭・ポスター発表の別）
Large-scale finite element simulationof earthquake ground motion and its
numerical verification aiming forenhanced earthquake damage
estimation

国内・外発表した時期の別（学会等名）
Ryoichiro Agata, Kohei Fujita,
Tsuyoshi Ichimura, Takane
JpGU - AGU Joint Meeting 2020: VirtualHori

2020/7/12

国内

2020/7/12

国内

セッション「ハイパフォーマンスコンピューティングが拓く固体地球科学の未来」の Discussion Forum Session
縣亮一郎, 藤田航平, 市村強, 堀 JpGU-AGU Joint Meeting 2020,(DFS）niteDFS
高峰
Virtual Meeting, 12－16 July 2020
[STT54-06] 地震被害想定に向けた地震
波動場の大規模有限要素計算とその数値検証（招待講演で口頭発表）

2020/7/13

国内

Finite element solver for large-scale
three-dimensional soil liquefaction
analysis

日下部亮太, 市村強, 藤田航平,
堀宗朗, Wijerathne Lalith

JpGU-AGU Joint Meeting 2020

2020/7/13

国内

Fast and scalable finite-element
solvers for three-dimensional urban
earthquake simulation

藤田航平, 市村強, 山口拓真, 堀
JpGU-AGU Joint Meeting 2020
宗朗, Lalith Wijerathne

2020/7/13

国内

地震被害想定に向けた有限要素法と有限 縣亮一郎, 藤田航平, 市村強, 堀 令和 2 年度土木学会全国大会年次学術講演
差分法による地震動解析の比較
高峰, 横田崇会

2020/9/11

国内

大規模数値シミュレーションによる地震発生から地震動・地盤増幅評価までの統 堀高峰, 市村強, 藤田航平, 縣亮
日本地震学会 2020 年度秋季大会
一郎
合的予測システムの構築とその社会実装のための取組

2020/10/29

国内

Trans-dimension インバージョン手法に 富田史章, 飯沼卓史, 縣亮一郎,
よる沈み込み帯における固着すべり分布 堀高峰, 市村強, 藤田航平, 山口 日本地震学会 2020 年度秋季大会
拓真
推定

2020/10/31

国内

Fast Scalable Implicit Solver with
Ichimura, T., Fujita, K.,
Convergence of Physics-Based
Koyama, K., Kusakabe, R.,
Simulation and Data-Driven Learning: Minami, K., Inoue, H.,
Toward High-Fidelity Simulation with Nishizawa, S., Tsuji, M.,
Digital Twin City.（ポスターセッショ Nishiki, T., Hori, M.,
Maddegedara, L. and Ueda, M.
ン）

The International Conference for High
Performance Computing, Networking,Storage and Analysis（SC20）.,SC20
Research Posters Display.,
Vurtual Meeting, 19 November 2020

2020/11/19

国外

高性能計算物理シミュレーションと人工知能の融合の試み―地震シミュレーショ 市村強ンを例に―

2020 年度第 1 回計算科学フォーラム（HPCIC 計算科学フォーラム）オンライン開催､2020 年 11 月 24 日

2020/11/24

国内

地震応答解析の解析能力向上のための高 藤田航平, 市村強, 堀宗朗,
Lalith Maddegedara
性能計算アプローチの提案

⽇本地震⼯学会・⼤会－2020

2020/12/3

国内

2020/12/11

国内

High Performance Computing & AI に市村強
よる大規模地震シミュレーション

第 62 回地盤工学シンポジウム オンライン
2020/12/18
開催､2020 年 12 月 18 日

国内

データ駆動型高性能計算物理シミュレー
市村強ションへ向けて

スーパーコンピューティング技術産業応用協議会セミナー

国内

Development and application of large
Takane Hori, Tsuyoshi
scale numerical simulations for
Ichimura, Kohei Fujita,
earthquake generation, ground motion
Ryoichiro Agataand non-linear soil amplification

JpGU - AGU Joint Meeting 2020: VirtualComparison of ground motion analyses
Ryoichiro Agata, Kohei Fujita,
using finite element and finite
The 3rd International Conference on
Tsuyoshi Ichimura, Takane
difference method aimed for enhanced
COMPSAFE 2020
Hori, Takashi Yokota
earthquake damage estimation

89

2021/3/10

２．学会誌・雑誌等における論文掲載
掲載した論文（発表題目）

発表した場所

発表者氏名（学会誌・雑誌等名）

Low-Order Finite Element Solver with Yamaguchi, T., Fujita, K.,
Ichimura, T., Naruse, A., Wells, Proceedings of the Platform for Advanced
Small Matrix-Matrix Multiplication
J.C., Zimmer, C.J., Straatsma, Scientific Computing Conference
Accelerated by AI-Specific Hardware
for Crustal Deformation Computation T.P., Hori, M., Maddegedara, L.
dU d N
HP3C 2020: Proceedings of the 2020 4th
Fast Multi-Step Optimization with
Ichimura, T., Fujita, K.,International Conference on HighDeep Learning for Data-Centric
Yamaguchi, T., Hori, M.,
Performance Compilation, Computing andSupercomputing
Wijerathne, L. and Ueda, N.
Communications

発表した時
期

国内・外
の別

2020/6

国外

2020/6

国外

Data-Driven Approach to Inversion
Analysis of Three-Dimensional Inner
Soil Structure via Wave Propagation
Analysis

Yamaguchi, T., Ichimura, T.,
Fujita, K., Hori, M.,
Wijerathne, L. and Ueda, N.

Computational Science – ICCS 2020. ICCS
2020. Lecture Notes in Computer Science,
vol 12139. Springer, Cham.

2020/6

国外

Data Assimilation in Volcano
Deformation Using Fast Finite
Element Analysis with High Fidelity
Model

Murakami, S., Yamaguchi, T., Computational Science – ICCS 2020. ICCS
2020. Lecture Notes in Computer Science,
Fujita, K., Ichimura, T.,
Maddagedara, L. and Hori, M. vol 12139. Springer, Cham.

2020/6

国外

本プロジェクトで進めてきた AI と物理
シミュレーションを融合した都市地震解 市村強析が紹介された（取材協力）

ITmedia/TechFactory オンライン記事 令和
2 年 8 月 4 日付け

2020/8

国内

Development of element-by-element
kernel algorithms in unstructured
finite-element solvers for many-corewide-SIMD CPUs: Application toearthquake simulation

Fujita, K., Horikoshi, M.,
Ichimura, T., Meadows, L.,
Nakajima, K., Hori, M. and
Maddegedara, L.

Journal of Computational Science, Volume
45, 101174

2020/9

国外

Development of regional simulation of
Kusakabe, R., Fujita, K.,
seismic ground-motion and induced
Ichimura, T., Yamaguchi, T.,liquefaction enhanced by GPU
Hori, M. and Wijerathne, L.
computing

Earthquake Engineering and Structural
Dynamics (EESD)., Special Issue Papers:
Regional Seismic Ground-Motion
Simulation and Observation withEngineering Applications.
First published 10 November 2020

2020/11

国外

The International Conference for High
Performance Computing, Networking,Storage and Analysis（SC20).,11th
Ichimura, T., Fujita, K., Hori,
M., Maddegedara, L., Ueda, N. Workshop on Latest Advances in Scalable
Algorithms for Large-Scale Systems
and Kikuchi, Y.
(ScalA20)
Vurtual Platform, 12 November 2020

2020/11

国外

地震本部ニュース 令和 2 年（2020 年）冬号
令和 3 年 1 月 18 日発行 第 13 巻第 3 号

2021/1

国内

2021/2

国外

2021/2/9

国外

A Fast Scalable Iterative Implicit
Solver with Green's function-based
Neural Networks

災害軽減へ向けて理学・工学をつなぐ計
市村強
算科学 ― 東京大学地震研究所 ―High-fidelity nonlinear low-orderunstructured implicit finite-element
seismic simulation of important
structures by accelerated element-byelement method

Fujita, K., Koyama, K.,
Minami, K., Inoue, H.,
Journal of Computational Science,
Nishizawa, S., Tsuji, M.,
Volume 49, February 2021, 101277
Nishiki, T., Ichimura, T., Hori,
M. and Maddegedara, L.

High-fidelity elastic Green's functions Takane Hori, Ryoichiro Agata,for subduction zone models consistent Tsuyoshi Ichimura, Kohei
with the global standard geodetic
Fujita, Takuma Yamaguchi,
reference system
Takeshi Iinuma

90

EARTH PLANETS AND SPACE, 73, 41(2021). https://doi.org/10.1186/s40623-02101370-y