メールマガジン No.51（2022 年 3 月）

JACM 賞受賞者セミナーシリーズ開催報告
岡田 裕
第 1 回 JACM 賞受賞者セミナーの開催報告を致します．
2021 年 12 月 10 日（金）夕刻に第一回を開催しました．
これは，JACM 賞受賞者の皆様の最新の研究成果を JACM会員の皆様と共有する目的の下，JACM（Japan Associationfor Computational Mechanics）の新しい取り組みとして開始したセミナーシリーズ第 1 回目になります．今後，年 3 回程度，最近の JACM 賞受賞者の方々にご講演をお願いし，オンライン形式で開催していきます．なお，新型コロナウ
イルス感染症終息後は対面・オンラインハイブリッド形式
として、講演終了後に簡単な懇親会を開催することで，ばらんな意見交換と研究交流の場になることを期待しています．
第一回目は下記の内容にて開催いたしました (写真 1) ．
日時：12 月 10 日(金) 18:30-20:00
場所：オンライン開催
1．開会あいさつ JACM 会長 岡田裕（東京理科大学）
2．ご講演：2021 JACM Young Investigator Award 堀田大介
氏（気象庁気象研究所気象観測研究部 主任研究官）
ネスト可能でマルチグリッド法の利用に適した二次元球面上の格子系と求積法

3．2019 JACM Young Investigator Award 杉本振一郎氏（八戸工業大学 准教授）
階層型領域分割法および並列電磁界解析への適用
4．閉会あいさつ JACM 副会長 塩谷隆二（東洋大学）
第 1 回目は約 25 名の方々にご参加頂きました．私個人の感想ですが，堀田様の全球気象モデルに関するご講演内
容は，固体力学系を軸とした工学系研究者である私にはなじみのない内容でしたが，大変興味深く拝聴させて頂きました．地球温暖化とその影響で発生する気象変動が懸念される昨今，堀田様の全球気象モデルによる高解像度気象予測がますます重要になると感じました．また，杉本先生の人体モデルの高周波電磁界解析では，1,300 億自由度という，途轍もなく大規模な問題を，15 分で解析できるという，スーパーコンピュータと，それを活用する技術に驚かされました．
ご講演下さった堀田様と杉本先生に心より感謝申し上げます．ご講演内容の詳細は，本号に堀田様と杉本先生からご寄稿頂いています．そちらをご覧ください．
第 2 回は 2022 年 3 月 17 日 17 時より開催予定です．

写真 1 参加者集合写真- 1 -JACM 会長（東京理科大学）

ネスト可能でマルチグリッド法の利用に適した二次元球面上の格子系と求積法
堀田
筆者は 2021 年 7 月 29 日にオンライン開催された JACM総会にて JACM Young Investigator Award を受賞する名誉に預かり、2021 年 12 月 10 日には JACM セミナーにて受賞
対象研究の概要を紹介する機会を頂きました。本稿ではセミナーでの講演の要旨をご報告します。

大介（気象庁気象研究所）

いのですが、スペクトル法を用いる気象モデルでは格子空間から波数空間への球面調和変換における求積法の制約からガウス緯度と呼ばれる緯度方向に不規則な格子配置を採用することが普通です。著者らの目標を達成するにはまず、厳密な球面調和変換を可能としつつ入れ子状の構造をもった球面上の格子を考え出す必要があります。

1. はじめに気象庁をはじめ世界の気象センターが発表する天気予報は、気象現象を大気の流体現象として捉え、その時間発展を流体力学の方程式を数値的に積分することで作成されています。このような流体力学の数値計算に基づく気象予報を数値予報と呼んでいます。数値予報に用いられる気象モデルを流体ソルバーとしてみた場合、以下のような特徴があります：(1) 球面のジオメトリー（素朴な緯度経度格子では極が特異点になる）(2) 「硬さ」stiffness（解に高速な音波・外部重力波が含まれるが、これらの波は気象現象に影響がないので解く必要がない）(1)の境界条件は、球面上の変数を球面調和級数に展開するスペクトル法により自然に表現できます。また(2)の硬さ問題は semi-implicit 法で解決することができますが、semi-implicit 法から帰結する楕円問題(ヘルムホルツ方程式)を解くときも、波数空間ではラプラシアンが対角化されることから球面調和展開によるスペクトル法が適しています。
こうした事情から全球気象モデルでは球面調和変換によるスペクトル法が広く使われていますが、スーパーコンピュータの大規模並列化が進んだ現在、そのスケーラビリティの確保が将来にわたって可能かが懸念されています。
2021 年現在、数値予報では数百ノードに渡る数千〜１万程度の CPU コアを使うことが一般的です。この規模の計算ではノード間の通信のコストが高くなりますが、球面調和変換の実行には全対全の通信が必要になるため、モデルの解像度が高くなるといずれ通信が律速となることが予想されます。このスケーラビリティ問題を解決するためには、スペクトル法への依存の度合いを減らしていくことが必要です。しかしながら、気象庁の全球気象モデルは気象庁業務の基幹モデルとして下流に多くのユーザを抱えており、急激に仕様をガラッと変えることは簡単ではありません。また気候変動予測（積分時間数十〜数百年）用の低解像度から中期天気予報（積分時間１週間程度）用の高解像度まで、様々な水平解像度を同じモデルでサポートする必要があるため、高解像度での通信コスト削減だけに特化した離散化手法を採用することも望ましくありません。
そこで著者らは、スペクトル法に基づく現行のモデルを、マルチグリッド法を実装可能となるよう拡張し、スペクトル法とマルチグリッド法のハイブリッド・モデルとすることを目指すこととしました。これは semi-implicit 法を用いる気象モデルの計算で最も厄介な楕円問題を解く際、マルチグリッド法が並列化性能と収束の速さの面で有効なためです。マルチグリッド法では粗い格子と細かい格子を行ったり来たりするため、内挿に伴う変換誤差や計算コストを避けるため格子が入れ子状になっていることが望まし
- 2 -

2.

新しい格子系と求積法の提案
厳密な球面調和変換が可能な入れ子状の球面上の格子は、これまで著者の知る限り知られていませんでした。そこで著者らは球面調和変換に Clenshaw-Curtis の求積法を用いる新しい格子系を考案しました。球面調和変換で一般的に用いられる Gauss の求積法では緯度方向の積分において被積分関数を sin(緯度)のルジャンドル多項式に展開し、項別に解析的に積分しますが、Clenshaw-Curtis の求積法では sin(緯度)のチェビシェフ多項式に展開します。こうすることで分点が緯度方向に等間隔に並ぶため、緯度方向に入れ子状で、かつ sin(緯度)の多項式を数値的に厳密に積分できる格子を構成できます。
この緯度方向の格子配置を用いて、球面を 24 個の三角形に分割した格子を構成すると（図 1）、各緯度円上の経度方向の格子間隔も等間隔にすることができ、東西方向には高速フーリエ変換、南北方向には Clenshaw-Curtis の求積法を適用することで厳密な球面調和変換を可能としつつ、粗い格子（図 1 左）の格子点の緯度・経度方向の中点に新たに格子点を配置することで細かい格子（図 1 右）を生成できる、入れ子状の格子系を構成することができます。
著者らはこのような格子系を提案し、24 面体 Clenshaw‐Curtis 格子と名付けました。詳細は割愛しますが、この格子系には、極付近に格子点が集中しない、各セルの面積が一様に近い（特に、太線で区切られた最も粗い格子の 24個の三角形は面積が等しい）、といった、スペクトル法のみならず差分法や有限体積法にとっても有利な性質があります。

図1: 提案する24面体Clenshaw‐Curtis格子。球面を太線で示された境界により、互いに面積の等しい24個の三角形に分割し、それぞれの三角形をさらに分割する。粗い格子(左)の三角形のセルを4つの三角形にすることで細かい格子(右)を構成でき、この操作を繰り返すことで任意の解像度の格子を生成可能。
3.

求積法の性質の数値的な確認
提案する手法と現行の手法の重要な違いは南北方向の求積法の違いのみですが、気象モデルや気候モデルでは長い時間積分の間に数千回・数万回もスペクトル変換を繰り返すため、求積法の性質の僅かな違いにも注意が必要です。
そのため、提案した手法を気象モデルに組み込む前に、提案する南北方向の求積法の性質、特にエイリアシング（切断波数を超える波数を持つ成分が低波数成分に混入する現象）に伴う誤差の特性を入念に確認しました。その結果、南北方向の分点数を切断波数の２倍以上取ることで旧来のガウスの求積法と遜色ない精度を確保できることを確認しました。
4.

気象モデルを用いた実験
求積法や球面調和変換の数値的な性質の違いを入念に確認した上で、2 次元及び 3 次元のモデルに実際に提案手法を実装し、標準的なテストケースの設定を用いて現行手法との違いを確認しました。図 2 に代表例を示すとおり、提案手法と現行手法の結果の差は非常に小さく、見た目では違い確認することができません。規格化した相対的な差は 10-4 程度以下のオーダーでした。

可能としつつ格子点法ベースのマルチグリッド法とも親和性の高い格子系を提案しました。また提案した格子系を実際に気象モデルに導入し、スペクトル法による時間積分で現行のモデルと同等の結果が得られることを確認しました。
本研究の最終的な目標は、「今回提案した格子上で格子点法とスペクトル法を併用可能とし、楕円問題を解く際に粗い格子ではスペクトル法、細かい格子では格子ベースの繰り返し法を用いる格子・スペクトル・ハイブリッドのマルチグリッド法を用いることにより、任意の解像度、任意のノード数で高速に動作する気象モデルを実現すること、しかも現行の気象モデルから枠組みを大きく変え図に滑らかに移行すること」です。ここで紹介した成果はこの壮大な計画の最初の一歩に過ぎません。本研究の成果を足がかりに、最終目標の実現に向け着実に開発を進めていきたいと考えています。
謝辞
今回受賞の対象となった研究[1]の共著者である気象庁
数値予報開発センターの氏家将志さん、気象研究所の吉村裕正さん、京都大学の石岡圭一教授には多くのご助言をいただきました。また、国立環境研究所の八代尚さんには推薦の労を取って頂きました。心より感謝申し上げます。
本研究のような気象学分野の研究を、工学分野の研究者が多い JACM で評価して頂いたことを大変光栄に思います。気象学分野は理学部出身の研究者が多く、気象モデルの開発者は工学分野の最新の研究を十分にフォローできていないのが現状かと思います。今回頂いたご縁を大切に、気象・気候分野での数値計算法の研究者と他分野の計算力学のコミュニティーとの交流を深められるよう努めていきたいと考えています。

図2: 温帯低気圧の発達を模した3次元モデルのテストケースにおける提案手法（緑の等値線）と現行手法（黒の等値線で低解像度、シェードで高解像度の結果を示す）の積分結果の比較。積分開始9日目の850hPa面に沿った気温を示す。３つの結果はほぼ一致し、等値線は重なっていて黒線は見えない。

5.

まとめと今後の展望
気象モデルの大規模並列環境下におけるスケーラビリティ問題への解決策として、本研究では、球面調和変換を参考文献[1] Hotta and Ujiie (2018, QJRMS) DOI: 10.1002/qj.3282

階層型領域分割法および並列電磁界解析への適用
杉本振一郎（八戸工業大学）
1. はじめに
12 月 10 日(金)にオンラインで開催された第 1 回 JACM賞受賞者セミナーにて「階層型領域分割法および並列電磁界解析への適用」のタイトルで講演させていただいた．この講演では，まず階層型領域分割法の概要についてお話しし，続いてその並列電磁界解析への適用についてお話しした．本稿では，これらについて概説する．
2. 階層型領域分割法
2.1. インターフェース問題
一般的な領域分割法において領域を分割する目的は，分散メモリ並列のためにデータを分割することである．そのため解くべき問題(有限要素方程式 等)は変化しない．一方，階層型領域分割法では元の問題の自由度を，領域分割によって生じる領域間のインターフェース(IF)自由度(𝑢 ，図 1)に静的縮約した IF 問題(図 2)を考える．まず IF 問題- 3 -を並列反復法で解き，次にそれぞれの小領域で𝑢 をDirichlet 条件として領域内自由度𝑢 を求めることで，全体の解を得る．
2.2. Domain-by-Domain
サイズが大きく，元の問題よりも密なシュア補元行列𝑆を陽に作成しないことも階層型領域分割法の特徴である．
並列反復法に現れる𝑆の行列ベクトル積は，小領域でのローカルシュア補元行列𝑆 の行列ベクトル積の重ね合わせに帰着される(図 3)．さらに𝑆 の行列ベクトル積は小領域での有限要素計算の結果を用いることで得られる．これは小領域をスーパーエレメントと考えた Element-by-Elementであり，Domain-by-Domain と呼ばれている．
2.3. 階層型の領域分割
階層型領域分割法という名称は，計算時間やメモリ使用量の観点から最適な領域サイズを取れるよう，階層型のデータ構造を取ることに由来する．最適な領域サイズは電磁界解析では領域あたり 100～200 要素と比較的小さく，一般にノードやコアなどの演算ユニットの数に対して領域数は十分に多くなる．そこで全体を MPI プロセスと同数の part に分割した後， part 内で複数の subdomain に分割(図 4)した階層型のデータ構造を取る．これにより，subdomain を最適な領域サイズにすることができる．この階層型の領域分割では，part 内における subdomain に関する演算の反復処理を共有メモリ並列にすることで，容易にハイブリッド並列を実現できることもメリットの一つである．

part 0part 1sub-sub-domain 0

domain 0

sub-sub-domain 1domain 1

図 4. 階層型の領域分割．
3. 並列電磁界解析
3.1. 電磁界解析(有限要素法)の特徴
電磁界解析で対象となるのは Maxwell 方程式である．構造や熱伝導と同様に，実対称行列を解く定常問題，非定常問題の他に，複素対称行列を解く準定常問題が存在する．
商用電源や電磁波は本来非定常問題として扱うべき対象であるが，正弦波状に変化することから時間微分項𝜕⁄𝜕𝑡を 𝑖𝜔 (𝑖：虚数単位，𝜔：角周波数[rad/s])で置き換えることで複素行列を一度だけ解く準定常問題に置き換えることができる．準定常問題では正弦波の腹と節に当たる時間での値が得られる．非定常問題として扱う場合に比べて大幅に計算時間を短縮できるが，ソルバには実・複素双方の行列を解く機能が求められる．
また 3 次元の電磁界解析では，スプリアス解の回避やゲージ問題の解決を目的として辺要素が用いられるのも特徴の一つである．構造や熱伝導で用いられる節点に自由度がある節点要素と異なり，辺要素では辺に自由度を持つ(図 5)．

domain 0Domain
Decomposition

domain 1: Interface DOFs(𝑢 ) : Inner DOFs(𝑢 ): DOFs(𝑢)図 1. 領域分割と自由度の分類．

𝑓

元の問題：𝐾𝑢𝐾
𝐾𝐾
𝐾

領域分割
𝑢𝑓𝑢𝑓

静的縮約
IF 問題：𝑆𝑢𝑔𝑆𝑆

𝑅 𝑆 𝑅 , ：シュア補元行列
𝐾𝐾

𝐾𝐾

図 2. 自由度の静的縮約．
𝑹は自由度を制限する 0-1 行列．

DOFs𝑦

図 5. 節点要素(左)と辺要素(右)．

𝑆𝑥
𝑅 𝑆 𝑅 ∙𝑥𝑦𝑆𝑥

𝐾𝐾 𝑥

𝐾 𝑧

𝑅 𝑆𝑥
𝐾𝐾𝐾

一般に，電磁界解析の有限要素方程式は悪条件である．
例えば静磁場問題では，Coulomb ゲージを考慮した弱形式は正定値ではない(図 6 (1))．またひどく悪条件であるため反復法が収束せず，しかも対角成分に 0 があるためピボッティング付きの直接法で解かなければならない．一方，
Lagrange の未定乗数𝑝は 0 なので，𝑝を無視することで条件数が改善し，反復法で解くことのできる弱形式(図 6 (2))
を得ることができる．しかし不定性があるため直接法では解けない．階層型領域分割法では IF 問題を反復法で解くため，静磁場問題に階層型領域分割法を適用する際には図6 (2)の弱形式を用いなければならない．また subdomain の有限要素計算も直接法より精度が悪い反復法で解かなければならないため，IF 問題の収束性悪化という問題を抱えることになる．

𝑅 𝑦,𝑥

𝐾 𝑧,

𝐾 𝑥：小領域での有限要素計算.
図 3. シュア補元行列の行列ベクトル積．
(Domain-by-Domain)- 4 -3.3. 大規模電磁界解析例[3][4]高周波誘電加温による癌の温熱療法の効果を定量的に評価するために情報通信研究機構が公開している数値人体モデル(図 7)の高周波電磁界解析に取り組んでいる．温熱療法の対象となる癌の直径は数 mm 程度であるため，0.1 mm 以下の解像度(1 兆自由度以上)の電磁界解析を可能にすることを目指している． 2016 年に解像度 0.5 mm (160億自由度)の解析を，2019 年には解像度 0.25 mm (1,300 億自由度)の解析を行った．前者は東京大学情報基盤センタ
ーの Oakleaf-FX 4,800 ノードで約 10 分，後者は OakforestPACS 8,192 ノードで約 15 分で求解できた．現在は大規模可視化(図 8)に取り組みつつさらに電磁界解析の解像度を向上させるとともに，熱伝導との連成解析を進めている．

(1) Coulomb ゲージを考慮した弱形式
𝜈 𝑟𝑜𝑡 𝐴, 𝑟𝑜𝑡 𝐴∗
𝑔𝑟𝑎𝑑 𝑝, 𝐴∗
𝐽, 𝐴∗ ,0.
𝐴, 𝑔𝑟𝑎𝑑 𝑝∗・正定値ではない・悪条件：直接法では解ける・𝑝は 0 であることがわかっている(2) Lagrange の未定乗数を消去した弱形式
𝜈 𝑟𝑜𝑡 𝐴, 𝑟𝑜𝑡 𝐴∗
𝐽, 𝐴∗ .
・条件数改善：反復法で解ける・不定性がある：直接法では解けない
図 6. 静磁場問題の弱形式とその特徴．
𝐴：磁気ベクトルポテンシャル[Wb/m]，
𝑝：Lagrange の未定乗数，𝜈：磁気抵抗率[m/H]，
𝐽：電流密度[A/m2]，上付きの∗：任意の試験関数．
3.2. 階層型領域分割法の特徴を活かした取り組み
3.2.1. subdomain への直接法の適用[1]前述のとおり不定性のある方程式を用いた場合にはsubdomain の有限要素計算に直接法を用いることができず，IF 問題の収束性が悪い．そこで IF 問題には反復法で解くことのできる図 6 (2)の弱形式を用いつつ，subdomain の有限要素計算には直接法で解くことのできる図 6 (1)の弱形式を用いることで，IF 問題の収束性を大幅に改善した．弱形式の違いによる齟齬は，𝑝が 0 であることから，IF 上の𝑝に関する自由度に 0 の Dirichlet 条件を与えることで解消している．
3.2.2. 移動体を含む対象の高速解析[2]モータや発電機など，内部に移動する部品を持った回転機も電磁界解析の対象である．一般に有限要素法で移動体を含む対象を扱うのは困難であるが，移動体が移動しても，移動しない固定体との接合面上の要素面が一致するメッシュを用いることで回転機の解析が行われている．しかしこのことが回転機の並列解析を困難にしており，1 桁程度しか解析時間を短縮できない要因となっていた．そこで固定体と移動体を個別に領域分割し，接合面上の自由度を IF自由度にすることで解析時間を大幅に短縮した．この手法
では階層型領域分割法のアルゴリズム本体には一切変更を加えず，接合面上の自由度に関する通信テーブルを時間進展に伴って交換するだけで済む．

- 5 -

図 7. 数値人体モデルの外観．

図 8. 右上半身の磁場の等値面．
(解像度 0.5 mm，160 億自由度)
4. おわりに
本稿では第 1 回 JACM 賞受賞者セミナーにて講演させていただいた内容を概説した．この講演が聴講された皆様の今後の活動の一助になったならば幸いである．
参考文献[1] 金山 他, IEEJ B, 129-8, pp.1018-1024, 2009.
[2] Sugimoto, et al., JASSE, 4-1, pp.99-116, 2018.
[3] Sugimoto, et al., MEL, 3, p.16-0067, 10p, 2017.
[4] Sugimoto, et al., JSST2019, pp.290-295, 2019.

2022 JACM Awards 候補者の募集
岡田 裕

日本計算力学連合（JACM）は，計算力学分野における顕著な功績および業績をあげた研究者を表彰する 3 種類の JACM Awards の候補者を募集します．JACM 会員におかれましては，候補者を自薦他薦で奮ってご推薦下さい．
受賞者には，2022 年 8 月開催予定の 2022 JACM 総会において表彰予定です．
***************************************************The JACM Computational Mechanics Award
日本計算力学賞（3 名以内）
計算力学の広い分野での顕著な研究業績，ソフトウエア開発，計算技術開発を行った研究者に対して与えられる．
The JACM Young Investigator Award
日本計算力学奨励賞（3 名以内）
計算力学分野で顕著な業績及び研究を行った 40 才以下（表彰年 4 月 1 日現在で 40 才未満）の研究者に与えられる．

JACM 会長（東京理科大学）

推薦書に記載して頂く項目は以下の通りです．
1. 推薦しようとしている Award の名称
2. 候補者の氏名，所属・住所，e-mail アドレス
（奨励賞候補者は生年月日も記載のこと）
3. 推薦者の氏名，所属・住所，e-mail アドレス
4. 主な受賞歴を含む経歴（最大 10 行以内）
完全なリストである必要はありません．最近のもの，あるいは最も重要なポストを記載してください．
5. 候補者の最も主要な功績あるいは業績の簡潔な記述（500 字以内）．特に，その Award の候補者として推薦する理由がわかるように記載してください．
推薦状のフォーマット兼例文は，以下の URL にあります．
https://ja-cm.org/Japanese/Award/index.html
推薦書は，2022 年 3 月 31 日（木）までに e-mail にて，次のアドレスにお送りください．
送付先：award@ja-cm.org

The JACM Fellows Award
日本計算力学連合フェロー賞（5 名以内）
計算力学分野で顕著な業績を上げ，JACM へのサポート，および IACM 関連国際学会に貢献した研究者に対して与えられる．

なお，候補者は主に日本国内において活動した研究者（外国人も含む），候補者及び推薦者は，応募時点でJACM 会員であることが必要です．

過去の受賞者は，下記 URL で一覧できます．
https://ja-cm.org/Japanese/Award/past.html**************************************************
編集責任者
塩谷- 6 -隆二（東洋大学）