ISSN 2436-0066

Vol.55 No.1

2022-1኱つᶍ⛉Ꮫィ⟬ࢩࢫࢸ࣒㛵㐃᱌ෆ㸺኱つᶍ⛉Ꮫィ⟬ࢩࢫࢸ࣒㛵㐃ᴗົࡣࠊࢧ࢖ࣂ࣮ࢧ࢖࢚ࣥࢫࢭࣥࢱ࣮ᮏ㤋ෆࡢ᝟ሗ㒊᝟ሗᇶ┙ㄢࡀᢸᙜࡋ࡚࠸ࡲࡍࠋ㸼KWWSVZZZVVFFWRKRNXDFMS㝵  䝃䞊䝡䝇᫬ 㛫ᖹ᪥ィ⟬ᶵ฼⏝඲⯡䛻㛵䛩䜛┦ㄯ 䡚኱ุ䝥䝸䞁䝍䚸฼⏝⪅➃ᮎ➼䛾 䡚฼⏝฼⏝⪅ㄯヰᐊ

 ⮬㈍ᶵ

䡚

ᒎ♧ᐊ  ศᩓ䡶䢙䢇䢛䡩䡬䡼༤≀㤋 㐃⤡䛟䛰䛥䛔䚹Ṕ௦䛾኱ᆺィ⟬ᶵ➼䛾ᒎ♧䠖䡚䠖⥲ົ䛻㛵䛩䜛䛣䛸䡚ಀ䞉ᐊྡ฼⏝┦ㄯᐊ୍㝵㟁ヰ␒ྕ䠄ෆ⥺䠅 
HPDLO  
VRGDQ#FFWRKRNXDFMS┦ㄯဨ୙ᅾ᫬

ぢᏛᕼᮃ䛾᪉䛿ඹྠ฼⏝ᨭ᥼ಀ䜎䛷䛤୺䛺䝃䞊䝡䝇ෆᐜ⥲ົಀ  
FFVRP#JUSWRKRNXDFMS఍ィಀ  ఍ィ 䛻㛵䛩䜛䛣䛸䚸㈇ ᢸ㔠 䛾
FFNDLNHL䠜䡃USWRKRNXDFMS ㄳồ䛻㛵䛩䜛䛣䛸䡚

䠄䠅୕ ඹྠ฼⏝ᨭ᥼ಀ  㝵ཷ㻌 ௜ XNHWXNH#FFWRKRNXDFMS฼⏝ᡭ⥆䛝䚸฼⏝┦ㄯ䚸ㅮ⩦఍䚸䝷䜲䝤䝷䝸䚸ぢᏛ䚸䜰䝥䝸䜿䞊䝅䝵䞁䛻㛵䛩䜛䛣䛸䡚

ඹྠ◊✲ᨭ᥼ಀ

  
UVVHF#FFWRKRNXDFMSඹྠ◊✲䚸ィ⟬ᶵ䝅䝇䝔䝮䛻㛵䛩䜛䛣䛸

䡚

䝛䝑䝖䝽䞊䜽ಀ  䝛䝑䝖䝽䞊䜽䛻㛵䛩䜛䛣䛸
LQHWZRUN#JUSWRKRNXDFMSᅄ◊✲㛤Ⓨ㒊
㝵  ஬㝵 ➃ᮎᶵᐊ

 

3& ➃ᮎᶵ䠄; ➃ᮎ䠅䡚

䡚㸦 㸧ෆࡣᮾ໭኱Ꮫෆࡢࡳࡢෆ⥺㟁ヰ␒ྕ࡛ࡍࠋ㟷ⴥᒣ࣭ᕝෆᆅ༊௨እ࠿ࡽࡣ㢌࡟  ࢆຍ࠼ࡲࡍࠋᮏㄅ䛾ྡ๓䛂㻿㻱㻺㻭㻯䛃䛾⏤᮶㻌᫛࿴  ᖺ䛻ᮾ໭ᆅ༊䛾᭱ึ䛾㟁Ꮚィ⟬ᶵ䛸䛧䛶䚸ᮾ໭኱Ꮫ㟁Ẽ㏻ಙ◊✲ᡤ䛻䛚䛔䛶᏶ᡂ䛥䜜䛯䝟䝷䝯䝖䝻䞁ᘧィ⟬ᶵ䛾ྡ๓䛷 6(1$& 6(1GDL$XWRPDWLF&RPSXWHU 䛛䜙䛸䛳䛶࿨ྡ䛥䜜䛯䚹㻌SENAC Vol. 55, No. 1（2022. 1）

― 1 ―

［巻頭言］

サイバーサイエンスセンターの改組について
菅沼拓夫

東北大学サイバーサイエンスセンター長
令和 4 年 1 月 1 日付で、東北大学サイバーサイエンスセンターの研究開発部の構成変更
による改組を行いましたので、ご報告させていただきます。
近年、国内外において、研究データを中心に据えてその解析結果に基づき研究プロセスを進める「データ駆動型研究」が、大学における研究方法論の重要な柱として位置づけられてきております。例えば内閣府の「第 6 期科学技術・イノベーション基本計画」(2021 年 3 月)においては、その具体的な実現方策として、研究データの管理・利活用促進、スマートラボ・
AI 等を活用した研究の加速と、そのための研究施設・設備・機器の整備・共用、研究 DX を進めることによる新たな研究コミュニティ・環境の醸成が必須であるとされています。これまで本センターは、長年に渡り、大規模科学計算システムの利用と研究をサポートする体制を整え、また、その役割を高度化するために全国共同利用・共同研究拠点としてサービスを広く提供しつつ、その高度化に関する研究開発を推進することで、大規模計算科学の多くの研究者の皆様のコミュニティをご支援させていただいてきました。大規模科学計算システ
ムは、データ駆動型研究においてもデータの解析を行うための心臓部に当たる重要な役割を担うことは言うまでもありませんが、その利用を支援する体制を維持しつつ、令和 4 年度から始まる第 4 期中期目標・中期計画期間に向けて、データ駆動型で研究を推進する皆様にも本センターのサービスをより安心して便利かつ高度にご利用いただけるようにすること
で、さらに幅の広い利用者の皆様方の参加によるコミュニティの拡大・活性化に繋げるべく、センター機能のさらなる充実を検討してきたところです。
以上の検討結果から、このたび新たに、データ駆動型研究の中核となる「データ」に着目し、研究ビッグデータの安全かつ高度な利活用に資する研究部として、データプラットフォーム研究部、情報セキュリティ研究部の 2 研究部を新設、合わせて、現状に即して既存研究部の役割・名称変更を伴うセンター改組を、令和 4 年 1 月 1 日付で実施いたしました(図 1)。
本センターは、東北大学大型計算機センターが 1969 年に設置されて以来、平成 13 年の情報シナジーセンターへの改組、平成 20 年のサイバーサイエンスセンターへの発展的改組を経て現在の形となっておりますが、それに次ぐ今回の大きな改組となります。
各研究部の役割に関しては、具体的には、データプラットフォーム研究部では、大規模研
究データに対応した高度データプラットフォーム、データレイク等のデータストレージ関連のハード、ソフトに関する研究開発と、研究ビッグデータの効果的な獲得・蓄積・有効利用のためのデータプラットフォームの企画・運用・管理支援を行います。また、情報セキュ
リティ研究部は、研究データ等の流通や研究データクラウドにおける高度情報セキュリテ
ィ技術に関する研究開発と、情報セキュリティ強化システムの企画・運用・管理支援の役割

― 2 ―

SENAC Vol. 55, No. （1 2022. 1）を担います。加えて、情報通信基盤研究部から名称を変更した情報サービス基盤研究部では、サービスレイヤに注力し、大学 DX における情報サービス基盤とその高度な利用に関する研究開発と、クラウドサービスを含めた情報サービス基盤および各種情報サービスの企画・運用・管理支援を推進します。さらに、先端情報技術研究部から名称を変更したサイバーフィ
ジカルシステム研究部は、アプリケーション分野としてサイバーフィジカルシステムに着目し、その医療・福祉・介護への応用に関する研究開発と、サイバーフィジカルシステムにおける情報の適正な利用に関する教育・啓蒙活動を進めます。今回はまずは組織の改革からスタートしておりますが、今後、新たなスタッフとともに、センターの関連設備の強化についても積極的に施策を検討していく所存です。
引き続き、新しい組織にて、情報基盤の整備と人材育成に取り組みつつ、本センターが我が国をリードする優れたサービスを提供し、学内外の関係者の方々、特に本センターを利用される皆様のご期待に応えられるよう全力で取り組んでまいりますので、今後ともご支援とご指導を賜りますよう、よろしくお願い申し上げます。

東北大学サイバーサイエンスセンター
研究開発部 (～令和3年12月31日)

東北大学サイバーサイエンスセンター
研究開発部 (令和4年1月1日～)

ネ

ネ

ッ

ト

ワ

ー

ク

研

究

部

ッ

ト

ワ

ー

ク

研

究

部

スーパーコンピューティング研究部

スーパーコンピューティング研究部

情

報

通

信

基

盤

研

究

部

情 報 サ ー ビ ス 基 盤 研 究 部

先

端

情

報

技

術

研

究

部

サイバーフィジカルシステム研究部
情 報 セ キ ュ リ テ ィ 研 究 部
デ ー タ プ ラ ッ ト フ ォ ー ム 研 究 部

高性能計算技術開発（NEC）共同研究部門

高性能計算技術開発（NEC）共同研究部門

クラウドサービス基盤（CSI）研究室

クラウドサービス基盤（CSI）研究室

図 1 令和 4 年 1 月 1 日付センター改組(研究開発部の組織改編)

SENAC Vol. 55, No. 1（2022. 1）― 3 ―[共同研究成果]［共同研究成果］

ナノアンテナの数値解析のためのブロックモーメント法の高速化
今野 佳祐 † , 陳 強 †東北大学大学院工学研究科 通信工学専攻 †

まえがき

1

近年，筆者らはナノアンテナの数値解析のためのブロックモーメント法の高速化を行ってきた．昨年度の報告では，ブロックモーメント法を用いてインピーダンス行列を数値計算する部分のベクトル化を実施し，その一部を高速化することに成功した [1]．その一方で，インピーダンス行列の数値計算を行う部分のうち，指数積分の数値計算のベクトル化が困難であったため，インピーダンス行列の数値計算アルゴリズム全体のベクトル化には成功していなかった．指数積分の数値計算では，精度のよい結果が得られるアルゴリズム

が引数の値や積分路によって異なる．そのため，引数の値や積分路によっては必ずしもベ
クトル化しやすい数値計算アルゴリズムを選択できないことがあり，ベクトル化には工夫が必要である．

本報告では，昨年度に引き続きブロックモーメント法のベクトル化を行い，インピーダンス行列の数値計算アルゴリズム全体のベクトル化を行った結果を示す．ベクトル化が困
難な指数積分の数値計算の部分は，精度を損なわない範囲でベクトル化に有利なアルゴリ
ズムに変更し，数値計算の高速化を図った．数値計算例を示し，高速化アルゴリズムの有効性を明らかにする．

体積積分方程式を用いたブロックモーメント法

22.1

線状モノポールセグメント間の自己・相互インピーダンス

さて, 図 1 の観測セグメント上の電流分布を表す基底関数として, 以下の区分正弦関数を利用する [2]，[3]．

sinh jk0 (ea − ea1 ).
sinh jk0 |ea2 − ea1 |(1)sinh jk0 (sa2 − sa ).
fs (s ) =sinh jk0 |sa2 − sa1 |(2)

fe (ea ) =同様に，波源セグメント上の電流分布を表す基底関数として，以下の区分正弦関数を用いる．

a― 4 ―

SENAC Vol. 55, No. （1 2022. 1）

次に, (1) 式を連続の式に代入することによって, 観測セグメント上の線電荷密度を以下のように得る.

∇ · fe (ea ) + jωρe (ea ) = δ(ea − ea2 ){}1sinh jk0 (ea − ea1 )
a
a
a
−jk0
+ δ(e − e2 ) .
ρe (e ) =jω
sinh jk0 |ea2 − ea1 |(3)

同様に, (2) 式を連続の式に代入することにより, 波源セグメント上の線電荷密度を以下のように得る.

∇ · fs (sa ) + jωρs (sa ) = −δ(sa − sa1 ){}1sinh jk0 (sa2 − sa )
a
a
a
jk0
− δ(s − s1 ) .
ρs (s ) =jω
sinh jk0 |sa2 − sa1 |(4)

ここで, 以下に示すセグメント間の自己・相互インピーダンスの一般的な表示式に (1)-(4)

式を代入して整理すると, 端電荷を持つ線状モノポールセグメント間の自己・相互インピー

ダンス表示式が, 単積分の形で (5) 式のように得られる [4].
∫ ea2 ∫ sa2∫ a∫ a−jk0 Rijjωµ0 e2 s2jω
e−jk0 Rij a a
a
a e
a
a
Zij =
fe (e ) · fs (e )ds de +ρe (ea )ρs (sa )ds de
4π ea1 sa1Rij
4πε0 ea1 sa1Rij= C A I A + CB I B + CC I C − C D I D .

(5)√ただし, (5) 式中の CA ∼ CD , IA ∼ ID , そして Rij はそれぞれ，Rij = (ea − sa )2 + (sb2 ) + (sc2 )2 ,
2
2 ∑
∑
a
a
e−jk0 (ne1 +ms2 ) mnIpq , CB = j sin k 15sa −sa , IB =CA = sin k ea −ea7.5sin k sa −sa , IA =
0| 2
0| 2
0| 2
1|
1|
1|
p=1 q=12
2∑∑aa ∫ v(sa ,ea ) −z
IBp , IBp = −mejk0 (me2 −ms2 ) v(sa2,ea2) e z dz, CC = j sin k 15ea −ea , IC =ICq , ICq =
1 2
0| 2
1|
p=1
q=1
−jk0 R21aa ∫ w(sa ,ea ) −z
30
, ID = e R21である．なお，IA の表示式に含−nejk0 (−ne1 +ns1 ) w(sa1,ea2) e z dz, CD = jk
0
まれる Ipq は以下のようになる．
∫ a∫ aaa
1 + mn e2 s2 ejk0 (ms +ne −Rij ) a a
Ipq = jk0
ds de
mn
Rij
ea
sa
1
1
∫ v(sa1 ,eah ) −z
∫ w(sah ,ea2 ) −z2∑eea (m+n)h
jk0 ea(n+m)
jks
0dz − e hdz=(−1) {−e ha ,ea )
a ,ea )
z
z
v(s
w(s
1 h
h 1
h=1
∫ x(sah ,ea2 ) −z
∫ y(sah ,ea2 ) −zee
+ ejk0 ju0
dz + e−jk0 ju0
dz}.
a
a
z
z
x(sa
y(sa
h ,e1 )
h ,e1 )
1

1

(6)

ただし，v(sai , eah ) = jk0 (−msai +meah +Rhi ), w(sah , eai ) = jk0 (−neai +msah +Rih ), x(sah , eai ) =√, d=
jk0 (−msah −neai +Rih +ju0 ), y(sah , eai ) = jk0 (−msah −neai +Rih −ju0 ), u0 = d 1+mn
1−mn√
(sb2 )2 + (sc2 )2 , m = (−1)p−1 , n = (−1)q−1 であり，(6) 式において, 1 + mn = 0 なら

Ipq = 0, 1 − mn = 0 なら最後の 2 項を無視する. Ipq ，IB ，IC の表示式において, z を含む全ての積分は 2 つの指数積分の差を取ることで数値的に積分できる．なお，モノポールセ

グメントの配置が垂直の場合も同様の定式化が必要だが, そちらは省略する．

Ipq ，IB ，IC の表示式において，指数積分の数値計算が最も時間のかかる部分である．

ナノアンテナの数値解析のためのブロックモーメント法の高速化

― 5 ―

ea
i

d

( s2b )2  ( s2c ) 2 (e2a , 0,0)


b
b
c
c


a
s
s

,
s
s
(e1 , 0, 0)  1
2 1
2


eb

O
d

j

( s2a , s2b , s2c )( s1a , s1b , s1c )ec
図 1: 端電荷を持つ一対の線状モノポールセグメントの一例と座標系.

2.2

ブロックモノポールセグメント間の自己・相互インピーダンス
2 wec

ea

le i

O

ec


( s2b )2  ( s2c ) 2 (e2a , 0,
0)  dbccs1b s
,ss(e1a , 0, 0) 
2 1
2

2 web
eb
c
2ws
( s2a , s2b , s2c )dls j( s1a , s1b , s1c )2 wsb

図 2: 端電荷を持つ一対のブロックモノポールセグメントの一例と座標系.
図 2 に, 端電荷を持つ一対の平行なブロックモノポールセグメントを示す. 観測・波源

モノポールセグメントにおける電流の向き, 端電荷の符号と位置は前項の図 1 と同じである. 観測ブロックモノポールセグメントの電流分布を表す基底関数は, 以下のように定義される.

fe (ea ) =1 sinh jk0 (ea − ea1 ).
4web wec sinh jk0 |ea2 − ea1 |(7)

1 sinh jk0 (sa2 − sa ).
4wsb wsc sinh jk0 |sa2 − sa1 |(8)

同様に, 波源モノポール上の電流分布を表す基底関数は以下の式で定義される.

fs (sa ) =これらの基底関数は, 先頭の係数部分を除けばそれぞれ (1) 式, (2) 式と一致する. 従って,

ブロックセグメント中央を通る線状モノポールセグメント間の相互インピーダンス Zij を(5) 式によって求め, Zij を 4 重積分すると, ブロックモノポールセグメント間の自己・相互インピーダンス表示式が以下のように得られる.
∫ web ∫ wec ∫ sb1 +wsb ∫ sc1 +wsc
1
block
Zij dsc dsb dec deb .
Zij =16web wec wsb wsc −web −wec sb1 −wsb sc1 −wsc(9)

Zij の計算には単積分が含まれるため, (9) 式は 5 重積分を含む表示式である. この 5 重積

分は，座標変換を用いて 3 重積分の和 (互いに垂直なブロックセグメントであれば 4 重積

― 6 ―

SENAC Vol. 55, No. （1 2022. 1）分) に変換することができ，数値計算の収束性および計算時間が改善する [5]．詳細は省略

するが，本報告ではこの 3 重積分の表現を用いて数値計算を行う．ナノアンテナの数値解
析では，これらの式を用いて誘電体ブロックセグメント間の自己・相互インピーダンス行
列要素を計算し，得られた行列方程式を解くことで電流を求める．

2.3

指数積分のベクトル化

誘電体ブロックセグメント間の自己・相互インピーダンスの数値計算において，指数積分の数値計算には最も長い計算時間がかかる．したがって，指数積分の数値計算のベク

トル化は自己・相互インピーダンスの数値計算の高速化を図る上で必須であるが，NECSX-Aurora TSUBASA 向け数学ライブラリ（NEC Numeric Library Collection）には複素数を引数とする指数積分のライブラリがない [6]．そこで，本報告では指数積分の数値

計算のベクトル化を独自に試みた．本報告での指数積分の数値計算に用いたサブルーチン
は，[7] からダウンロードしたものであるので，そのベクトル化を実施した．

指数積分はいわゆる複素積分である．被積分関数は原点に極を持ち，積分路がブランチ

カット (=負の実軸) を横切ると異なるリーマン面へと移動する多価関数である．したがって，積分路が極の近くを通ったり，ブランチカットを横断するときの積分には注意が必要である．このような指数積分を数値的に求める方法には，漸近展開や級数展開，数値積分の他，連分数展開を用いる方法が挙げられる．これらの中で，漸近展開や級数展開，数値積分を用いたアルゴリズムのベクトル化は容易であるが，連分数展開を用いたアルゴリ

ズムはループ内の変数に依存関係ができてしまうのでベクトル化が困難である．したがって，指数積分の数値積分の高速化だけを考えるのであれば，連分数展開以外のアルゴリズ
ムを用いればよいことになる．一方で，数値計算の精度を考えると，引数の値によって数値計算のアルゴリズムを変えるのが最もよいことになり，引数の値によっては連分数展開
を用いることは避けられない．

本報告では，以上の議論を踏まえ，引数の値の範囲と取り得る積分路によって指数積分の数値計算のアルゴリズムを変える．互いに平行な誘電体ブロックセグメント間の自己・相互インピーダンスの数値計算では，指数積分の引数は純虚数になり，積分路が原点を横断することがない．そこで，このときの指数積分の数値計算は全て数値積分で実施する．
一方，互いに垂直な誘電体ブロックセグメントの自己・相互インピーダンスの数値計算

では，指数積分の積分路が原点あるいはブランチカットを横断することがある．そこで，
このときは場合分けを行い，その結果に基づいて指数積分の数値計算アルゴリズムを変える．指数積分の引数が複素数のとき，被積分関数の特異性が強い可能性がある (要するに，積分路が極やブランチカット付近を通る可能性がある)．したがって，このときは精度を優先し，数値計算は連分数展開と級数展開を用いて行う．指数積分の引数が純虚数のとき，被積分関数の特異性が弱ければ (要するに，積分路が極を横断していなければ) 高ナノアンテナの数値解析のためのブロックモーメント法の高速化― 7 ―速化を優先し，その数値計算は数値積分で行う．反対に，被積分関数の特異性が強ければ精度を優先し，その数値計算は連分数展開と級数展開で行う．このように，引数の値の範囲と取り得る積分路によって指数積分の数値計算のアルゴリズムを変えることで，精度を損なうことなくベクトル演算率の向上と高速化を図る．なお，数値計算に級数展開や連分数展開を用いる箇所にもループアンローリングを施し，可能な限りベクトル演算率を向上させた．

3

数値計算結果100000

CPU time [sec.]

10000
1000
100
0.001N2

10
1

1

10

100
1000
Number of unknowns N

10000

図 3: 数値計算時間．

銀でできた 100 × 100 × 100 nm3 の立方体を垂直入射の TM 波で励振し，その電流を数値的に求めた．銀の複素比誘電率は Drude モデルによって与え，周波数は 100 THz とした [8]．ナノダイポールアンテナはブロックダイポールおよびブロックモノポールセグメ
ントに分割した．

このときのセグメント数 N に対するトータルの数値計算時間を図 3 に示す．数値計算

時間は概ね N 2 に比例しており，これはインピーダンス行列の数値計算時間のオーダーと一致する．一般的に，モーメント法の数値解析では行列方程式を解くための数値計算時間

が N 3 に比例して増加するので，トータルの数値計算時間は N 3 に比例して増加することが多い．一方，今回の数値計算では，インピーダンス行列の数値計算に多重積分が含まれており，そのせいでインピーダンス行列の数値計算がトータルの計算時間の大部分を占めている．その結果，トータルの数値計算時間が N 2 に比例して増加している．このよう― 8 ―

SENAC Vol. 55, No. （1 2022. 1）に，誘電体の数値解析では，自己・相互インピーダンスの数値計算の高速化が必須である
ことが分かる．

プログラム全体のベクトル演算率は 95.3%であり，指数積分の箇所以外のベクトル化を実施した昨年度のプログラムのベクトル演算率 2.7%と比較して大幅にベクトル演算率が向上した．自己・相互インピーダンス 1 つあたりの数値計算時間は 0.0016 sec. であり，昨年度のプログラムの数値計算結果の 0.43 sec. と比較して 200 倍以上の高速化を実現した．

4

まとめ

本報告では，ナノアンテナの数値解析のための体積積分方程式を用いたモーメント法の高速化を行い，指数積分の数値計算の箇所のベクトル化に成功した．その結果，プログラ

ム全体の高速化が図れ，数値計算の大幅な高速化に成功した．今後は，様々なモデルで本プログラムの数値計算精度を検討し，アンテナの設計にも応用していく予定である．

謝辞

本研究成果は, 東北大学サイバーサイエンスセンターのサポートを受けたものである．

スタッフをはじめ，関係各位に深謝する．

参考文献[1] 今野 佳祐，陳 強，“ナノアンテナの数値解析のためのブロックモーメント法の高速化，” 東北大学大規模科学計算システム広報 SENAC vol. 54, no. 1, pp.25-31, Jan.

2021.

[2] J.H. Richmond and N.H. Geary, “Mutual impedance of nonplanar-skew sinusoidal
dipoles,” IEEE Trans. Antennas Propag., vol.AP-23, no.3, pp.412-414, May 1975.
[3] 陳 強，ザイ フイチン，袁 巧微，澤谷 邦男，“誘電体に対するガラーキンモー
メント法 -端部電荷を考慮した直方体モノポール間の自己・相互インピーダンスの単積分化-，” 信学論 (B)，Vol. J91-B，No. 9，pp. 926-939，2008 年 9 月．

[4] M.A. Tilston and K.G. Balmain, “On the suppression of asymmetric artifacts arisingin an implementation of the thin-wire method of moments,” IEEE Trans. Antennas
Propag., vol.38, no.2, pp.281-285, Feb. 1990.
[5] A. Köksal and J.F. Kauﬀman, “Mutual impedance of parallel and perpendicularcoplanar surface monopoles,” IEEE Trans. Antennas Propag., vol.39, no.8, pp.12511256, Aug. 1991.

ナノアンテナの数値解析のためのブロックモーメント法の高速化― 9 ―[6] NEC Numeric Library Collection 2.3.0 ユーザーズガイド，https://www.hpc.nec/

documents/sdk/SDK_NLC/UsersGuide/asl/c/ja/index.html, 2021 年 8 月 16 日ア

クセス．

[7] Special Functions, “https://people.sc.fsu.edu/~jburkardt/f77_src/special_
functions/special_functions.html,” 2021 年 8 月 6 日アクセス．

[8] M. Nafari, and J. M. Jornet, “Modeling and performance analysis of metallic plasmonic nano-antennas for wireless optical communication in nanonetworks,” IEEE
Access, vol. 5, pp. 6389-6398, 2017.

― 10 ―

SENAC Vol. 55, No. （1 2022. 1）［共同研究成果］

細線化ならびに多層化されたプリント基板の FDTD 解析における高速化の検討
⻑野⼯業⾼等専⾨学校

春⽇

貴志

プリント基板は⾼速⼤容量のデータ処理や多機能化に伴い，多層化ならびに配線の細線化が⾏われている．有限差分時間領域（FDTD）法を⽤いた電磁界解析を実⾏する場合，プリント基板構造により最⼩の単位セルサイズが決定するが，セルサイズに対して基板⾯積が広いため解析が膨⼤となる．本研究では，MPI による計算の⾼速化について検討した．セル数によるベクトル化率についての検討では，256 の倍数の前後においてベクトル⻑は増えるが，演算時間とベクトル化率に⼤きな影響がないことがわかった．また，ノード数の⽐較では，ノード数を増やすことで概ね増えたノード数だけ⾼速化が実現できていることが分かった．SX-ACE と SX-Aurora TSUBASA (以下、SX-Aurora)の⽐較では，SX-ACE に対して SXAurora は約 4〜5 倍の⾼速化が実現できることが明らかとなった．

1．はじめに
5G や 4k・8k 放送などの普及により，⼤容量データを処理する電⼦機器の性能は向上している．電⼦機
器の⾼性能化に伴い，プリント基板内の信号伝送は⾼速化している．プリント基板内の信号伝送における信号品質（Signal Integrity: SI）の改善や，ノイズ放射（Electromagnetic Interference: EMI）を低減させる
ためには，設計時に回路シミュレータや電磁界シミュレータにより SI や EMI の検討[1〜4]を⾏った上で，設計へフィードバックする．しかしながら，⾼性能化や⾼機能化により，プリント基板の部品点数は増加
しており，信号配線が増加するために基板の多層化ならびに配線の細線化が進んでいる．現在，プリント基板の層間の距離や線路幅，線路間隔は約 0.1mm オーダーが多い．このような基板の電磁界解析を⾏う
ためには，構造的な⼨法を正しく解析モデルに反映させるため，線路構造よりも⼩さなセルでモデル化する必要がある．⼩さなセルサイズに対し，基板の⾯積は⼤きいため，電磁界解析を⾏う際には⼤規模解析となる．
本報告では，厚みに対して表⾯積が遙かに⼤きいプリント基板を対象とした電磁界解析において，MPIを⽤いた有限差分時間領域（FDTD）解析の⾼速化について検討した．また，今回は SX-ACE と SX-Auroraの更新時期に共同研究に採択されたこともあり，それぞれのコンピュータ同⼠の性能の⽐較も⾏った．
2．解析対象図 1 に解析対象のプリント基板モデルを⽰す．今回は，計算の⾼速化に関する内容であるため，プリン
ト基板は簡易的なモデルとした．FDTD 解析における単位セルサイズはx=y=z=1mm である．時間ス
テップは 10,000 ステップとした．

細線化ならびに多層化されたプリント基板のFDTD 解析における高速化の検討― 11 ―

図 1 プリント基板

16 ノード
図 2 MPI による解析イメージ

32 ノード

図 3 16 ノードと 32 ノードの区分

FDTD 解析を行う場合，MPI プロセス数の違いによる解析の高速化を行う．図 2 に MPI プロセスを用いた 8 ノードの空間分割のイメージ図を示す．FDTD 法の場合，プリント基板とその周辺を取り巻く自由空間，解析の最外殻に吸収境界条件（ABC）として Perfect Matched Layer（PML）を設ける必要がある．
プリント基板と自由空間の電界と磁界は，各ノード共に共通のアルゴリズムにより解析が行えるが，PMLは最外殻に適用されるため，ノード毎に PML の配置を変更する必要がある．このため，MPI プロセス数が増えるほど PML の解析プログラムの行数が増える．
比較対象として，16 ノードと 32 ノードのイメージを図 3 に示す．16 ノードは x 方向に 4 分割，y 方向に 2 分割した．32 ノードは x と y 方向をそれぞれ 4 分割した．z 方向に分割数を増やさない理由として
は，プリント基板の厚みは薄いため，z 方向のセル数が増加しない傾向があるためである．一方で，x とy 方向は基板の平面方向にあたり，基板面積が増えるほどセル数が増大する．

図 4 ベクトル化率の検討図

― 12 ―

SENAC Vol. 55, No. （1 2022. 1）

3．ベクトル化率の検討

250
98.8

4 の解析モデルにより解析を行った．x 方向のセル

Real Time [sec]

数は 2028 から 2068cell まで変化させた．x 方向を他

平均ベクトル長160の y や z よりも大きくすることで，x 方向をベクト

200

140

ル演算させる．演算には SX-Aurora を用いた．図 5

Real Time
平均ベクトル長に MPI 8 の解析における解析時間（Real Time）
，平
120

98.6

98.4

ベクトル演算率

均ベクトル長，ベクトル演算率の変化を示す．なお，

ベ ク ト ル 演 算 率 [%]

セル数に対するベクトル化率を検討するため，図

2040
2060
x方 向 セ ル 数

この値は Rank 0～7 プロセスの平均値である．ベク
トル演算率に注目すると 98% 以上になっており，

図 5 MPI 8 による演算時間とベクトル長，ベク

高効率な計算ができているといえる．またセル数がトル演算率（SX-Aurora）

2064 cell と 2068 cell のときを比較すると，2068 cell では平均ベクトル長が 28.8，ベクトル演算率が 0.12%
低下し，Real Time が 9.3 秒増加した．セル数が 2068 cell のときに Real Time が増加しベクトル演算率が低下したのは MPI 通信の処理によるものと考えられる．
図 6 に MPI プロセス数 16（以下 MPI 16）による演算時
250

解析空間の中央のノードには，吸収境界条件 PML の解析

100
Real Time
平均ベクトル長80

平均ベクトル長について，
x 方向セル数が 2048 cell
（=256

図6

2040
2060
x方 向 セ ル 数

MPI 16 による演算時間とベクトル98100なったために，9 個目のベクトル長の効率が落ちたためとときに Real Time が増加しベクトル演算率が低下したのは MPI 通信の処理によるものと考えられる．

Real Time [sec]cell ではベクトル演算率が低下した．セル数が 2044 cell の

80

200

して吸収層の厚みのセル数が少ないため，ベクトル長に影

97.6

97.4

Real Time
平均ベクトル長ベクトル演算率2040
2060
x方 向 セ ル 数

算率の結果を示す．MPI 16 と同じ傾向が得られた．MPI 8
る．PML は，電波吸収のために 9 層設けている．面に対

97.8

150
60

図 7 に MPI32 での演算時間とベクトル長，ベクトル演だけ異なる傾向になった理由としては PML の影響であベ ク ト ル 演 算 率 [%]250

考えられる．
また，ベクトル演算率については，x 方向セル数が 2044

100

長，ベクトル演算率（SX-Aurora）

cell では平均ベクトル長が低下した．これは 1 プロセス当たりの解析空間の x 方向セル数が 256 の倍数より大きく

97.8

ベクトル演算率28.1 秒増加した．

98

150

向セル数が 2044 cell になったとき，急激に Real Time が×8）までは高い平均ベクトル長を維持しているが，2052

平均ベクトル長については，PML の違いによる差は見られなかった．x 方

98.2

200

平均ベクトル長

ル長や演算率については，ノードにより異なる．演算時間

Real Time [sec]が，最外殻のノードと比較して少ない．このため，ベクト

120

98.4

ベ ク ト ル 演 算 率 [%]間と平均ベクトル長，ベクトル演算率を示す．MPI16 では，

図7

100

MPI 32 による演算時間とベクトル

長，ベクトル演算率（SX-Aurora）

細線化ならびに多層化されたプリント基板のFDTD 解析における高速化の検討― 13 ―

響が出たと考えられる．
4．ノード数による比較
図 8 に MPI のノ

3000

ード数を変化させ

MPI 8

ACE ， (b) が SXAurora による演算
時 間 で あ る ． SX-

2000

MPI 161000

はメモリ量が足り

0400
MPI 16200

MPI 32

MPI 32

Aurora の MPI8 の1400cell について

MPI 8

Real Time [sec]示 す ． (a) が SX-

Real Time [sec]た際の演算時間を600

ないため計算が実
行できなかった．

800

1000
1200
y方 向 セ ル 数

0

1400

800

1000
1200
y方 向 セ ル 数

1400(b)SX-Aurora(a)SX-ACE

図 8 ノード数による演算時間の比較結果

SX-ACE の 1400cell で比較すると，MPI8 に対して MPI16

50
48

では 43%，MPI32 では 25%に時間が短縮し，高速化が実て MPI16 では 60%，MPI32 では 35%に時間が短縮した．
ノード数が 2 倍，4 倍に増えたにも関わらず，演算時間が50%，25%にならなかった理由としては，分割処理やプロ
セス間通信によるオーバーヘッドの影響と考えられる．
SX-ACE と SX-Aurora との比較では，MPI 8 の 1000cellのとき，演算時間は SX-ACE が 2275.7 秒，SX-Aurora が426.4 秒であり，SX-ACE よりも 18.7% に時間が短縮した．同じく，MPI16 では 25.6%，MPI32 では 26.0%まで時間が短縮した．新しい SX-Aurora では，約 4～5 倍の高速

使用メモリ量 [GB]現できる．一方，SX-Aurora の 1200cell では，MPI8 に対し化が実現できた．

600

40

30

20
500

1000
y方向セル数1500

図 9 使用メモリ量の推定（SX-Aurora）

SX-Aurora の MPI8 でメモリ制限により計算が実行できなかったことについて考察する．SX-Aurora のVE のメモリ上限は 48GB である．図 9 に，MPI 8，16，32 の解析で使用したメモリ量を示す．解析規模に比例して使用するメモリ量が増加していく様子が読み取れる．1 つの VE に搭載されているメモリは48 GB であるため，これ以上メモリを消費するような解析はできない．MPI 8，16，32 で解析できる最大の空間の規模を予想すると，図 9 のメモリ使用量の予測の点線と横軸との交点を読み取り，x 方向を4200，z 方向を 400 と一定にしたとき
MPI 8 :

4200 × 1300 × 400cell

MPI 16 : 4200 × 1550 × 400cell
MPI 32 : 4200 × 1700 × 400cell程度と推定した．細線化した基板などの解析のため，これ以上のセル数のモデルを扱う場合には，使用す― 14 ―

SENAC Vol. 55, No. （1 2022. 1）

るノード数を増やして計算を分散させる必要があると考えられる．
5. まとめ
本報告では，厚みが薄く，面積が広いプリント基板構造のような偏った解析対象を想定し，SX-ACE な
らびに SX-Aurora により計算の高速化について検討した．その結果，セル数によるベクトル化率についての検討では，256 の倍数の前後においてベクトル長は増えるが，演算時間とベクトル化率に大きな影響がないことがわかった．また，ノード数の比較では，ノード数を増やすことで概ね増えたノード数だけ高速化が実現できていることが分かった．SX-ACE と SX-Aurora の比較では，SX-ACE に対して SX-Aurora
は約 4～5 倍の高速化が実現できている．ただし，SX-Aurora の VE のメモリ上限が 48GB であるため，解析規模が大きくなると実行できないことが分かり，大規模解析には MPI の数を増やすことが不可欠である．
今後，プリント線路の細線化や高密度化が進むと予想され，解析に求められる規模が 10 年前から飛躍的に大規模化している．このため MPI64 や MPI128 についても検討が必要である．
謝辞：
本研究は，令和 2 年度萌芽型課題「MPI によるプリント基板の電磁界解析法の高速化」として行った．高速化に対して，東北大学サイバーサイエンスセンターからの支援に感謝する．
参考文献[1] T. Kitazawa, R. Kitahara, T. Yamagiwa, J. Chakarothai, Y. Hayashi, T. Kasuga, "Basic Study on a Novel FDTD

Method Implemented Frequency Dispersion of PCB", 2021 Joint IEEE International Symposium on EMC+SIPI
and EMC Europe, TP-WED-8B-4, (2021).
[2] 北澤太基，春日貴志，林優一他，“FDTD解析を用いたFR-4基板の異方性が与える信号伝送への影響の検討”，電子情報通信学会環境電磁工学研究会，EMCJ2021-17，（2021.7）[3] 北澤太基，春日貴志他，“コネクタにおける信号伝送評価のための3D解析モデルの構築”，電子情報通信学会環境電磁工学研究会，EMCJ2020-67，(2021.1)．
[4] 北澤太基，春日貴志他，“周波数分散性を組み込んだFDTD法によるFR-4基板の信号伝送解析に関する基礎検討”，電子情報通信学会環境電磁工学研究会，EMCJ2020-44，(2020.10)．

SENAC Vol. 55, No. 1（2022. 1）― 15 ―［共同研究成果］

流体解析における時間発展計算の AI サロゲートモデルの追求――― 直感による究極的な単純化はどこまで有効か？
松岡 浩
技術士事務所 AI コンピューティングラボ
筆者らは、東北大学サイバーサイエンスセンターの共同研究公募制度により、平成 30 年度から
「リカレントニューラルネットワークによる高解像度流体解析コードの開発」を行ってきた。すでに、SENAC Vol. 54 No.1(2021.1)において、筆者がめざしている「“格子ガス法流体解析モデル”「流体計算をと“ニューラルネットワーク”の融合[[1]」について、その趣旨を述べた。すなわち、行う際に設定する格子点に“多層パーセプトロン”が内蔵されていると考え、格子ガス法に基づ
く流体計算の時間発展計算を微小整数だけを用いるニューロン計算で行い、その結果を疎視化す
ると、ナビエ・ストークス(NS)方程式を解いた結果とほぼ同等の結果を得る」という手法の可能性について述べた。この手法は、NS 方程式による時間発展計算をニューラルネットワークの計算で代替するという意味で、本稿では、“AI サロゲートモデル”と呼ぶことにする。
これからの“ものづくりをめぐる計算機シミュレーション”の究極的な目標は、マルチスケー
ル、マルチフィジクス、リアルワールド(不完全な計測情報でも機能するリアルタイムデータ同化
等)、カーボンニュートラル(温室効果ガス排出削減のための超低消費電力計算等)の同時実現であると思われる。“AI サロゲートモデル”は、この目標に近づく有力な１候補であると考えられる
が、まだまだ力不足である。ここでは、この“AI サロゲートモデル”を極力単純化することによって、計算機の設計や性能に対する要求条件を大きく緩和し、少しでも上記目標の達成を容易にする可能性を追求した。
なお、ここで述べる内容は、学術論文的なものの要点解説ではない。まだ、研究途上であるた
め、これまでの共同研究を通じて、筆者が考えるに至った“AI サロゲートモデル”の方向性を、読み物的に紹介するものとした。また、計算法の詳細は、過去の SENAC に掲載した内容[[1-5]とほ【補足】に記
とんどが共通であり、本稿では、AOBA-A で行った計算の代表的な計算条件を末尾の【[6] ”載するにとどめた。計算法の基礎理論は、Teixeira が考案した“４次元 FCHC54 速度モデル[に基づく“格子ガス法[[7-10]”である。筆者なりの理論解説は、文献[[11]に掲載している。
１．ものづくりをめぐる計算機シミュレーションへの要求の激化
ものづくりをめぐる計算機シミュレーションへの要求は、時代とともにハードルが急速に高まっている。
計算機性能の飛躍的な向上とあいまって、はじめに強調されたのは、「実験や試作を最小限の数にできる“高精度なマルチスケール and/or マルチフィジクスな事象解析”の実現」であった。
次に、Society 5.0 に代表される高度な IoT 時代の到来とともに“リアルワールド問題”への対応を可能にする“不完全な計測情報でも機能するリアルタイムデータ同化(計測データ等のシミュ
レーション計算への実時間反映)”も重要な課題になった。しかし、はじめの時点では、スーパーコンピュータによる大規模計算が必要な“高精度なマルチスケール and/or マルチフィジクスな事
象解析”と、計測情報を入力しながら超高速計算にフィードバックする“リアルタイムデータ同化”の同時実現は、技術的にハードルが高すぎる遠い目標であり、それぞれ別の達成目標とみなされていた。その後、特に東日本大震災を契機にして、津波等に対する防災対策強化などの観点から、大規模計算とリアルタイムデータ同化を両立させる方向の研究開発が盛んになったように思う。また、人工知能(主に、ニューラルネットワーク)の第３期ブームが到来し、深層学習(DeepLearning)等の多方面への応用が注目されるようになり、大規模計算の超高速化にも成果をあげつ― 16 ―

SENAC Vol. 55, No. （1 2022. 1）

つある。また、これらの動きに並行して、スーパーコンピュータの計算性能は、半導体の微細加工技術による性能向上が限界に近付きつつあるとは言え、着実に高速化を続けてきた。さらに、現時点では応用分野が限られているとはいえ、量子コンピュータの実運用も開始された。こうした背景から、すでに、“マルチスケール”、“マルチフィジクス”、“リアルワールド(不完全な計測情報でも機能するリアルタイムデータ同化等)”の同時実現は、単なる妄想ではなく現実味のある目標になってきている。
そして今、世界は、“カーボンニュートラル最優先の時代”を迎えた。ものづくりについては、製品の設計、試作から、その利用・廃棄・リサイクルに至るライフサイクル全体での温室効果ガ
スの排出量削減が要求される。ものづくりをめぐる設計・試作・製造・利用(運転)などのために活用する計算機シミュレーションに対しても、温室効果ガスの排出につながる火力発電の利用を抑制するため、“超低消費電力化”の要求が今まで以上に重視されることになろう。
以上は一般論であったが、ここで具体的な例をひとつとりあげてみる。
たまたま、電気新聞を見ていたところ、「2021 年 11 月 24 日、スペインのシーメンス・ガメサ・
リニューアブル・エナジーが、フィンランド最大の陸上風力発電所向けに、風車 69 基を受注した。」との記事があった。これは、ロータ直径が約 170ｍ、タワーの高さはおよそ 150ｍということである。風車のブレード各部における空力計算に必要な解像度を大きめに見積もって仮に cm オーダーであったとしても、ロータ直径とのスケールの違いは４桁ある。また、ウインドファーム
においては、ある風車の後流は、下流側の風車の存在やその風車の後流と重なり、さらに下流の風車に対しても複雑な後流となって到達する。精度の高い発電量の予測には、地形を考慮した風況予測に加えて、このような後流の状況をウインドファーム全体について解析する必要がある。
このため、少なくとも数 km～数 10km 以上のスケールまでカバーするマルチスケール解析が要求される。また、複雑に変化する後流の中に存在する風車の構造解析では、弾性変形をともなう強度解析や、寿命の長い安全運転を実現するための疲労解析を行う必要がある。すなわち、流体と複雑な構造物の間の連成を伴うマルチフィジクス解析が必要になる。さらに、いろいろな風況の出現頻度の予測ができたとしても、リアルワールドで生じる時々刻々の風況予測は不可能である。
結局、運転中における疲労損傷の蓄積が小さくなるような運転制御を実現するには、現地センタ
ーの計測情報と連携したリアルタイムシミュレーションが重要になってくる。そして、最後に、
カーボンニュートラルの実現のために導入されるウインドファームについて、その最適配置、強度設計、運転制御等に必要なコンピュータシミュレーションで大電力を消費してしまうのでは意
味がない。
すなわち、
“ものづくりをめぐる計算機シミュレーション”の今後の究極的な目標は、マルチス
ケール、マルチフィジクス、リアルワールド、カーボンニュートラルの同時実現である。
２．“AI サロゲートモデル”という方向性への期待“マルチスケール”、“マルチフィジクス”、“リアルワールド(リアルタイムデータ同化)”、“カ
ーボンニュートラル(超低消費電力)”という技術的なハードルは、いずれも極めて高い。スーパ
ーコンピュータの性能が著しく向上してきたとは言え、これらを同時に解決できるほどの実用的な方策は見つかっていない。現時点では、ハードウェアとソフトウェアに関する個々の技術を積み上げて、少しずつでも目標に近づく努力をするほかはない。言い換えれば、目標達成には、まだかなりの時間を要すると思われる。他方、それでも、“ものづくりをめぐる次世代の数値シミュ
レーション技術”として、何か抜本的に大きな効果のある対応策はないものだろうか？と思いをめぐらすのが人情である。
以上のような状況の中、“AI サロゲートモデル”という方向性は、かなり効果のある対応策であると考えられる。

流体解析における時間発展計算の
AI サロゲートモデルの追求― 17 ―

一般に、“AI サロゲートモデル”という言葉からは、どのようなものを想像されるであろうか？
ただし、ここでは、冒頭で述べたように、本稿では、「NS 方程式による時間発展計算をニューラル
ネットワークの計算で代替するモデル」という狭い意味で、この言葉を用いている。
いまどきの AI は、ほとんどが“深層学習(DL: Deep Learning)”のイメージであろう。従って、
通常想像されるイメージは、例えば、以下のようなものではないだろうか。
【通常想像される“AI サロゲートモデル”のイメージ例】
マルチスケール解析を克服するため、まずは、非常に大規模な計算ができるというスーパー
コンピュータの“Capability Computing 能力”を活かし、NS 方程式による高解像度な時間発展計算を行う。これにより、１回の計算でも、「時刻ｔにおける詳細な流速空間分布」と「時刻ｔ＋⊿ｔにおける詳細な流速空間分布」のセットを非常に多数得ることができる。
次に、中規模な計算であれば非常に多数の計算を並列実行できるというスーパーコンピュー
タの“Capacity Computing 能力”を活かして、上記の計算で得られた各時刻の流速空間分布(速度ベクトルの各成分を表す３つの数値の巨大な３次元配列)を深層学習(DL)により次元圧縮する。次元圧縮によって、必要な記憶容量を大きく削減できる。ただし、次元圧縮で得られた数
値配列は、ニューラルネットワークの計算によって機械的に学習されたものであるので、その数値分布をみても物理的な意味は理解できない。
最後に、次元圧縮した時刻ｔと時刻ｔ＋⊿ｔにおける流速空間分布のセットを入出力データとしてニューラルネットワークに学習させる。このときの入出力データセットは非常に多数あるが、次元圧縮されているので、ひとつひとつの学習に要する計算規模は比較的小さい。
以上の結果、次元圧縮された流速空間分布の時間発展計算法が、ニューラルネットワークの計算法として学習されたことになる。これが、NS 方程式による時間発展計算を代替する“AI サロゲートモデル”のひとつのイメージである。
なお、この“AI サロゲートモデル”を用いて時間発展計算を行えば、直接得られるデータは、
やはり意味のわからない次元圧縮されたデータである。しかし、時間発展計算の過程で、我々がそのスナップショットを知りたい時はいつでも、次元圧縮の逆計算によってもとの高解像度レベルで物理的な意味のわかる流速空間分布を復元することができる。
以上のような技術は、深層学習(DL)の登場によってはじめて実用性が認識された AI 技術である。今後とも、その様々な活用によってこれまでにない重要な成果を期待できるであろう。
しかしながら、“マルチスケール”、“マルチフィジクス”、“リアルワールド(リアルタイムデー
タ同化)”、“カーボンニュートラル(超低消費電力)”の同時達成は、極めて高いハードルである。
この克服を目指すのであれば、何かさらに抜本的なアイデア、特に、上記の例で示した通常の“AIサロゲートモデル”のイメージをはるかに超えた“究極的な単純化”が必要だと思われる。
そこで、筆者は、以下、“AI サロゲートモデル”を究極的に単純化した姿について考察した。
３．“AI サロゲートモデル”の究極的単純化
今回の考察では、以下の３つの究極的な単純化について検討した。
（１）異なる物理現象の表現を単純化できないか？⇒“サブマクロ粒子”の挙動で代替
ものづくりをめぐる現実の流体解析の問題では、純粋な流体解析だけが要求される場合はほとんどない。構造物との熱流体力学的な相互作用をはじめ、相変化や、燃焼などの化学反応も同時並列して計算するべき問題も多い。これらの解析では、異なった物理現象が相互に関係しあって時間発展していくものと認識され、それぞれの物理現象に対して成立する支配方程式や経験式を

― 18 ―

SENAC Vol. 55, No. （1 2022. 1）

同時に満足するように連成計算がなされる。
他方、「すべての物質は、原子・分子からできており、異なるように見える物理現象も、結局は、“ミクロ粒子”の挙動で説明できるはず。」という自然科学的な知見がある。これは、マルチフィジクス解析をシングルフィジクス解析で処理するという単純化の方向を連想させる。ただし、「原子・分子スケールの計算から積み上げて、ものづくりをめぐる通常の工学スケールでの解析結果を得るという方法は、計算すべきマルチスケールの幅を極端に高めてしまう。」という重大な欠点がある。このため、ここでは、“ミクロ粒子(＝原子・分子)”ではなく、仮想的な“サブマクロ粒子”の挙動による物理現象の代替を考える。このことが、“単純化されたサロゲートモデル”を追求するうえで、一番重要な最初の選択になると考える。
【注】ものづくりをめぐる流体解析を行う場合、最終的に求めたい解像度のスケールにおける速度分布は、通常、“マクロな流体塊”の重心速度分布である。例えば、格子ガス法の場合、仮想粒子がもつ運動量などの物理量を、疎視化するために設定したある時空間領域で合計(平均)して時空間のスケールを一段上げることにより、“マクロな流体塊”がもつ物理量を求める。ここでは、このように、「一段スケールアップ(疎視化)しただけでマクロな流体挙動が得られるようなスケールで仮想する粒子」という意味を強調するため、“サブマクロ粒子”と呼ぶことにした。
また、“サブマクロ粒子モデル”の採用は、マルチフィジクスをシングルフィジクス化できるというメリット以外に、以下のようなメリットがある。すなわち、もともと流体解析だけというシングルフィジクスの場合を考えても、流体が激しく動きまわり、界面をもって砕け散るような事象については、通常の CFD では計算処理が困難である。このような事象に対しても、“サブマクロ粒子”のモデルは、一般的に対応しやすいとされている。
なお、現在の“サブマクロ粒子”のモデルを概観すると、一番注目されている手法は、「粒子法」であると筆者は思う。実際、従来手法では解析が困難な様々なマルチフィジクス事象について粒子法による解析事例が多数発表されている。また、「格子ボルツマン法」もかなりの成果を積み上げてきており、「格子ガス法」を呑み込んだ感じもある。以上は、筆者の個人的感想であるが、「粒子法」、「格子ボルツマン法」、「格子ガス法」は、「一定規模の原子や分子の粒子集団を仮想的な“サブマクロ粒子”に対応させてその挙動解析を行う」という点は共通であり、多様な
マルチフィジクスへ応用のアイデアは、相互に活用することができるであろう。
（２）サブマクロ粒子の状態記述を単純化できないか？⇒とりうる状態の数を微小整数に抑制“サブマクロ粒子”は、もともとが“多数の原子や分子の集団”である。このため、それがもつ物理量としては、いろいろなものを想定できる。具体的には、サブマクロ粒子内におけるミクロ粒子の粒子数密度、質量、運動量、角運動量、エネルギー、電荷などの分布やその合計値などである。また、それぞれの物理量がとりうる数値は、一般には、巨大な整数値や実数値になるので、とりうる状態の数は、ほとんど無限である。
しかし、ここでは、究極的な単純化を行う場合を考察しているので、まずは、サブマクロ粒子の挙動を特徴づける物理量を、“運動量”と“エネルギー”に限定してみる。
【注 1】仮想粒子に“質量”を担わせないモデルは、マルチフィジクスへの応用可能性を拡げる
可能性がある。例えば、弾性体内部の応力については、質量を伴わない“運動量の流れ(テンソル)”と解釈できる。また、運動量を担う同じ仮想粒子が、流体と弾性体の境界で相互に運動量を交換できるので、流体と構造物の間のより自然なシミュレーションを実行できると思われる。

流体解析における時間発展計算の
AI サロゲートモデルの追求― 19 ―【注 2】仮想粒子に“質量”を担わせないモデルは、仮想粒子が存在しない状態「０」について、新しい解釈を可能にする。すなわち、「０」は、仮想粒子が存在しない状態を直ちに意味するのではなく、単に、仮想粒子が静止している状態も含む。この静止状態の物理的な解釈は、「多数のミクロ粒子の集団が滑らかに移動できずに微小振動しながら渋滞している状態」のイメージになる。
さて、前述した３つ手法(粒子法、格子ボルツマン法、格子ガス法)は、“サブマクロ粒子”の挙動計算方法について、それぞれ異なる特徴をもつ。ここでは、究極的な単純化という意味で、浮動小数点表示による実数演算が不要であるというメリットを重視し、“格子ガス法[[7-10]”に基づく“サブマクロ粒子”のモデルを選択した。
第５章に後述する試計算では、Teixeira が考案した“４次元 FCHC54 速度モデル[[6]”を使用している。このモデルでは、54 種類の“サブマクロ粒子(仮想粒子)”を設定する。その内訳は、エネルギーが０で静止状態の仮想粒子が６種類、エネルギーが１で運動量の向きが異なる仮想粒子が 24 種類、エネルギーが２で運動量の向きが異なる仮想粒子が 24 種類である。また、「ひとつの格子点には、同じ種類のサブマクロ粒子は、高々ひとつまでしか存在できない」という格子ガス法に特有な通常の仮定をする。この仮定によって、時間発展計算時の演算が著しく単純化され、記憶容量も極めて小さく抑えられる。
ここで重要なことは、「格子ガス法」の選択によって、サブマクロ粒子がとりうる状態の数を「54」という微小な整数に抑えるという単純化が可能になることである。
【注】「54」という数値は、通常の格子ガス法の応用としては、数値が大きいように思われるかもしれない。しかし、“４次元 FCHC54 速度モデル”は、格子ガス法が正しく機能するための条件(ガリレイ変換に対する不変性、応力テンソルの等方性、非圧縮性流体に関する NS 方程式への漸近性などの条件)を自然な形で満足できることが、Teixeira の論文[[6]によって示されている。この点を重視して、筆者は、専ら 54 速度モデルを利用している。
（３）近傍粒子の状態考慮を単純化できないか？⇒考慮範囲を最小近傍に限定し方向ごとに処理
オリジナルな格子ガス法では、各格子点にいろいろな向きから到着した仮想粒子は、格子点上で瞬時に衝突散乱を起こし、いろいろな向きに出発していく。このとき、到着粒子と出発粒子の間で、粒子数、運動量およびエネルギーの総和が保存される。すなわち、出発粒子が出ていく向きのパターンは、到着粒子の向きのパターンのみが考慮された衝突散乱規則によって決定される。
ただし、粒子数や運動量の保存則を満たす出発粒子のパターンが複数存在した場合には、その選択は等確率でなされる。
他方、格子ガス法の幅広い応用に向けて、相分離[[3]や相変化などを伴う特別なシミュレーショ
ンを行う場合がある。この応用では、最隣接格子点までか、場合によっては、格子点間距離の数倍遠方までの範囲で、仮想粒子間に働く力を考慮して出発粒子のパターンを決定するモデルが採用されている。なお、粒子法の場合は、“影響半径”の概念がこれに相当し、それぞれの粒子から“影響半径”以内の距離に存在する粒子からの影響を考慮している。このように、近傍粒子の状態考慮は、幅広い応用の実現にあたっては不可欠な選択となっている。
そこで、ここでは、衝突散乱直後の出発粒子は、近傍に存在する仮想粒子の影響を受けて最終的にたどり着く格子点が変化しうるモデルとした。
[6]”では、4 次元面心超立方体(FCHC)格子を
Teixeira が考案した“4 次元 FCHC54 速度モデル[

― 20 ―

SENAC Vol. 55, No. （1 2022. 1）

用い、各格子点上に存在できる仮想粒子の数は、前述のとおり、最大 54 個までである。仮想粒子が持ちうる速度を４次元方向の成分で表現すると、静止粒子(0,0,0,0)が６個までと、運動粒子が、エネルギー１の粒子：(±1,±1,0,0),(±1,0,±1,0),(±1,0,0,±1),(0,±1,±1,0),(0,±1,0,±1),(0,0,±1,±1)で 24 個、エネルギー２の粒子：(±1,±1,±1,±1),(±2,0,0,0),(0,±2,0,0),(0,0,±2,0),(0,0,0,±2)で 24 個となる。従って、ある格子点から１回の時間ステップで到達することができる１番近い格子点は、４次元空間において、エネルギー１の粒子に関して 24
個、エネルギー２の粒子に関しても同じく 24 個である。これら 48 個の最隣接格子点の状態は、それぞれの最隣接格子点に存在しうる 54 個の仮想粒子の存否情報(存在するか否かの情報)で表現できる。従って、ある格子点の近傍状態を求めるには、48×54＝2592 個の仮想粒子の存否情報を扱う必要があり、かなり複雑になることが予想される。
ここでは、究極的な単純化を目指す観点から、考慮する近傍の仮想粒子の範囲を極力せばめてみた。具体的には、最隣接格子点より手前で、格子点間距離の中央位置までを考慮すべき“最小近傍”とした。格子ガス法では、静止粒子以外のすべての仮想粒子は、時刻ｔから時刻ｔ＋⊿ｔの間に格子点間を同期して移動する。ということは、時刻ｔ＋⊿ｔ/2 の時点では、すべての運動
粒子は格子点間の中央付近を移動中と考えることができる。いま、時刻ｔ－⊿ｔ/2 の時点において、上記の“最小近傍”の範囲に存在しうる仮想粒子を考える。それは、①注目している格子点を時刻ｔ－⊿ｔに出発して離れていく仮想粒子か、あるいは、②注目している格子点に時刻ｔに到着する予定で近づきつつある仮想粒子、のいずれかである。従って、存否情報を扱う必要がある仮想粒子の数は、48 個の向きがあっても各向きに２個ずつなので、合計で 48×2＝96 個に激減させることができる。
また、ここでは、さらなる単純化を行う。具体的には、96 個の近傍粒子の存否情報を複合して扱うのではなく、１方向(互いに反対の２つの向き)ごとに、その方向の“最小近傍”を移動中の仮想粒子だけを観察する。この場合、１方向の“最小近傍”には、最大でも４個の仮想粒子しか存在しない。これら４個の最小近傍粒子の存否パターンがある特定の条件を満足した場合に、注目している方向のうち、ある一方の向きに飛び出そうとしている仮想粒子があれば、その速度を反転させて反対向きに出発させるという操作を確率的に行う。ここでは、このモデルを“仮仮想粒子反転連行モデル”を呼ぶ。
４．“仮想粒子反転連行モデル”の概要
液体の密度は、固体の密度とそれほど変わらない。従って、液体中では、固体中と同様に、粒子はいたるところに密集して存在する。ある粒子の近傍には、その粒子を取り囲むように粒子が存在するので、その粒子の自由な行動は許されない。このイメージをヒントにして、以下のように考える。
個々の格子点における到着粒子の衝突散乱によって、たまたまある向きに飛び出そうとする出発粒子があったとしても、その前後には仮想粒子、すなわち、ミクロ粒子の巨大集団が存在する。
前方に存在する仮想粒子の速さが遅ければ、出発粒子は、前方の仮想粒子に追突し、跳ね返されて減速される。後方にそれより速い仮想粒子が存在すれば、今度は、後方の仮想粒子に追突されて加速される。結果的に、出発粒子は、その前後に存在する仮想粒子による挟まれ具合いに応じて速さを変化させつつ、その流れの中に“連行”されてしまうであろう。
【注 1】上記のような連行過程のイメージにおいては、一般的に、「“速度の向き”が変化する揺動」も考えられる。しかし、ここでは、前述のとおり、考慮する仮想粒子の運動方向を１方向に限定している。また、“速度の向き”は、時間平均あるいは集合平均では変化しないことから、個々の“連行”においても、出発粒子の速度が他の方向に変化する揺動はないものと仮定した。

流体解析における時間発展計算の
AI サロゲートモデルの追求― 21 ―

しかし、上述した「連行過程における個別粒子の“減速”や“加速”のイメージ」をそのままモデル化することはできない。なぜなら、１方向に限定した場合、運動粒子が持ちうる速さの選択
は、運動エネルギーが１または２の二者択一であり、個別粒子に微妙な“減速”や“加速”を行わせることはできないからである。そこで、“仮想粒子”が「もともとは粒子集団である」ことを思い出し、集団としての“速さ”の微妙な変化は、「出発粒子がもつ速度をそのまま保持しながら同じ向きに連行される確率」と、「出発粒子の速度が前後の仮想粒子との相互作用によって反転させられ、反対向きに同じ速さで連行される確率」の比率を用いて表現することにする。
【注 2】「仮想粒子の完全なる速度反転が生じる」という本モデルの仮定について、「多数の粒子が同時に反転するという事象が実際には起こりにくい」と思われることから、不適切な仮定ではないか？という疑問がわくかもしれない。しかし、一般に、サブマクロ粒子(仮想粒子)の挙動に関する仮定は、それを疎視化する時空間の範囲で合計(平均)したときに、マクロな流体塊の重心挙動を適切に表現できれば問題ないとされるべきである。よく用いられる“バウンズバック”の境界条件は、この代表的な事例であろう。
ただし、個々の“仮想粒子反転連行操作”は、仮想粒子の運動量を変化させる。特に何かの外力を模擬しようとしている場合でなければ、疎視化時空間の範囲で運動量を合計(平均)したときに、総和の保存が必要である。そこで、今回のモデルでは、ある向きへの連行と逆向きへの反転連行を、必ず交互に行うこととした。これより疎視化スケールでは、運動量保存則が成立している。
５ ． 試 計算に よ る 単純化 の 効 果に関 す る 考察（ １ ） １方向 ご と の最小 近 傍 仮想粒 子 の 挙動の 考 慮 だけで 何 が できる か ？
「サブマクロスケールにおける仮想粒子の挙動に対する単純な操作が、疎視化によって得られるマクロスケールの流体挙動において、どのような変化を発現させるか？」を定性的に観察するため、以下の試計算を行った。
はじめに、比較の基準を作るため、各格子点における衝突散乱の後、ある向きに飛びだした出発粒子が、「最小近傍粒子の影響を受けずにそのまま出発した場合」の円柱後流の様子を計算した。
図１に示す。
その結果を図ここでは、静止流体中に置かれていた無限長円柱に対して、時刻ゼロで、円柱軸に垂直な方向に
急に流体が動き出した場合の過渡変化を計算している。円柱後流の渦領域は、長い尾の形状になって一時的に円柱直径の約７倍の長さまで伸びたが、その後ゆっくり縮んで最終的には、約５倍の長さで安定した。
次に、「１方向ごとに最小近傍粒子の存否状態
を調べ、それに応じて、出発粒子の飛び出す向き
が前方か後方かの確率を意図的に偏らせた場合」に生じる円柱後流の様子を比較した。
図 １ ．近 傍 粒 子の影 響 を 考慮し な い 通常の 格 子 ガス 法 に よる円 柱 後 流の変 化
１方向の最小近傍粒子は、注目している格子点から見ると、“ある向き”と“その逆向き”の両

― 22 ―

SENAC Vol. 55, No. （1 2022. 1）

方に、それぞれ到着直前の粒子と出発直後の粒子が存在しうる。従って、最大４個の仮想粒子が存在しうる。最小近傍粒子の存否パターンは、1 方向あたり 16 パターンということになる。
このうち、４パターンは、格子点から見て、前後に存在する最小近傍粒子の配置が全く同じになる。この場合、最小近傍粒子の配置パターンを考慮しても、前後に飛び出す仮想粒子の出発確
率の間に意図的な偏りを生じさせる合理的な理由はない。そこで、この４パターンを除いた 12 パ
ターンすべてについて、その配置パターンが生じたときに、その直後の格子点における衝突散乱において「前方に飛び出す粒子がなく後方に飛び出そうとする粒子があれば、それを５％の確率で反転させる」という操作を行ってみた。
図２に示す。12 パターンのうち、３パターンでは、すぐにカルマン渦が発生した。
この結果を図また、長い尾が２パターン、短い尾が１パターン、非常に短い尾が１パターン、残りの５パター
ンでは、レイノルズ数が非常に小さい場合に生じる“円柱前後の対称流れ”が直ちに現れた。
これらの結果を分析すると、「１方向ごとに存在しうる４個の最小近傍粒子の合計の運動量が、ある向きに大きく、かつ、それに続けて同じ向きに仮想粒子を反転加速する確率が高いほど、す
なわち、周辺の粒子の流れに“連行される”と想像できる３パターンの場合に、円柱後流は、長い尾からカルマン渦に移行している」ということがわかる。これは、レイノルズ数が増加したときの変化である。試計算では、円柱の寸法は不変であり、入口のマクロ流速と各部の粒子数密度もほぼ一定に保たれているので、流体の動粘性が小さくなったものと解釈できる。従って、“１方向ごとの最小近傍仮想粒子の挙動を考慮する超単純なモデル”であっても、流体粘性を変化させる効果を発現でき、これを流体粘性の制御法として活用できる可能性がある。
【注】ここで、興味深い点がある。通常、格子ガス法で流体粘性を下げるためには、出発粒子を運動量勾配が正になる向きにできるだけ出発させるという方策がとられる[[13]。ところが、上記のモデルの場合、仮想粒子を反転させることによって出発確率を増す向きが、格子点前後における運動量の勾配がゼロ(下図の第４行第３列)あるいは負の向き(下図の第４行第５列)であっても、カルマン渦が発生し、流体粘性の低下が観察されている。
右図の第１行目は、前方の最小近傍に存在しうる到
着直前粒子と出発直後粒子の配置を示し、第１列目は、後方の最小近傍に存在する
到着直前粒子と出発直後粒子の配置を示している。
これら前後の配置パターンが同時に成立したときに、衝突散乱後の仮想粒子が出発する向きを、意図的に、前方を向く確率を増加させた場合の円柱後流の様子を示している。
なお、シミュレーションの計算体系等の詳細を【補足】
に示す。

前方
後方
カルマン
渦へ発展
円柱前後
対称流れ

長い尾

カルマン
渦へ発展

非常に
短い尾

長い尾

円柱前後
対称流れ

短い尾

円柱前後
対称流れ

円柱前後
対称流れ
カルマン
渦へ発展

円柱前後
対称流れ

図２ ． ⊿t / 2 過去 の 最 小近傍 仮 想 粒子の 配 置 に応じ て 衝 突散乱 直 後 の仮想 粒 子
を 前 方へ反 転 加 速する 場 合 の条件 と そ の結果 生 じ る円柱 後 流 の変化

流体解析における時間発展計算の
AI サロゲートモデルの追求― 23 ―

なお、上記の渦発生過程は、サブマクロスケールで以下のように解釈できる。
「各格子点の最小近傍には、合計 96 個の仮想粒子が存在しうる。このうち、48 個の到着粒子は、別の格子点から時々刻々やってくるので、これらを含めた 96 個の仮想粒子の重心速度は、激しく揺らいでいるであろう。従って、この１格子点ごとの重心速度(＝最大 96 個の仮想粒子小集団の重心速度)は、その格子点を含むマクロな疎視化領域(流体塊)における平均値、すなわち、“マクロな流速”を中心にしてそれなりに揺らいでいると考えられる。
１格子点ごとの重心速度が、たまたまある向きに揺らいだ時に、そのタイミングに合わせて続
けて同じ向きに仮想粒子が出発するような仮想粒子の反転連行操作を確率的にしかも前方と後方に向けて交互に行えば、この揺らぎは拡大を続け、やがて、疎視化時空間領域を超えたマクロスケールでの流体振動を発現するであろう。この流体振動は、より大きなマクロスケールの流体領域の重心を中心にいろいろな方向に発展していくので、そこに流体振動の位相差や流体領域の並進運動があれば、“目に見える渦”という形で視覚されるようになる。」
このように、「マクロな流体挙動は、サブマクロスケールで生じる仮想粒子の“連行”の強度、
すなわち、『格子点間隔スケールでの仮想粒子小集団の重心揺動が続けて同じ向きに生じる確率』に支配される。」ということができる。
45
また、渦発生過程とマクロな速度勾配の関係は、以下の通りに解釈できるであろう。
「一般に、マクロなスケールにおいて、速度勾配が一定の値に維持されている状態を考えてみる。マクロな流体塊の重心速度は、みな主流であるマクロな流速と同じ向きを向いている。このとき、マクロな46
流体塊の中をサブマクロなスケールで観察すると、各格子点の“最小近傍”に存在する仮想粒子小集団の重心速度は、主流と同じ向きではなくマクロな主
流ベクトルを中心にしていろいろな向きに揺らいでいると考えられる。この揺動状況を分析すると、速47
度勾配がある場所では、到着粒子の揺動主方向と出発粒子の揺動主方向は互いに乖離していることがわかる。この状態に、前述の“仮想粒子反転連行操作”を施してやると、到着粒子が来る頻度が少ない方向における“出発粒子の揺動”が拡大され、やがて、疎視化領域のスケールに達する。
48
このとき、主流速度の速い側から遅い側へ移動する流体塊は、主流に垂直でその流体塊の移動と同じ向きのマクロな流速速度成分をもつとともに、移動先の主流速度より常に少し速い主流方向のマクロな流速速度成分をもつであろう。逆に、主流速度の遅
49
い側から速い側へ移動する流体塊は、主流に垂直で
上記とは逆向きのマクロな流速速度成分をもつとともに、移動先の主流速度より常に遅い主流方向のマ
クロな流速速度成分をもつであろう。すなわち、主流に垂直なマクロな流体塊の往復振動は、その中心が主流に乗って下流に移動しながら楕円運動する流体の流れ、すなわち、渦の発生を意味している。そして、“仮想粒子反転連行操作”を行う確率を上げれば、速度勾配が生じる場所からの、このような渦の発生頻度を上げることができる。」

50

― 24 ―

試計算では、円柱の境界層における速度勾配からの渦発生が、“仮想粒子反転連行操作”の実行確率を上げることによって、より盛んになる様子が観測された。
図３に円柱境界層から発生して拡大していく渦の例(スナップショットの時刻ステップ 45～56)を示す。
なお、このような円柱から発生する渦が、このあと、ローラ渦や縦渦をうみ、複雑に重なった構造をもつ後流へと変化していくことが知られている[[12]。

SENAC Vol. 55, No. （1 2022. 1）

5152（ ２ ） どこま で 乱 流を模 擬 で きるか ？
乱流のように激しく変化する流れを計算機で安定的に計算していく観点からは、格子ガス法には重要なメリットがある。格子ガス法の時間発展計算は、ビット演算または微小整数演算のみで行うことができ、時間発展計算の個々の過程では、浮動少数点演算に起因する誤差の発生を避けることができる。このため、どんなに激しい流れについても、途中で計算不能になることはなく、必ず、何かの計算結果を安定的に求めることができる。
また、乱流状態は、大きなレイノルズ数の条件下で生じる。そして、多くの工学的なものづくり解析において、非常に大きなレイノルズ数条件下での解析が要求される。ここでは、どんな大きなレイノルズ数の解
析も理論上は実現可能なように、まず、流体粘性を負にできるシミュレーション方法を探求した。“仮想粒子反転連行モデル”の適用確率を調整することによって負の流体粘性を発現できることが、試計算によって定性的に確認できている。なお、格子ガス法による“負の粘性発現”を乱流解析に利用するアイデアについては、Rothman の論文[[13]を参照されたい。

53

54

55

56

図 ３ ． 円柱境 界 層 からの 渦 の 生成
さらに、実際にシミュレーション計算を実行する場合には、“負の流体粘性を発現させる方法”と“正の流体粘性を発現させる方法”を組み合わせて、模擬したい乱流状態に対応した小さい正の流体粘性を安定的に維持する具体的な方法が必要である。
図４に示す。
この実現を目指した試計算の例を図図２
具体的には、「衝突散乱直後の仮想粒子がその前後の流れにのみこまれる“連行参入”」が図の第２行第３列に示す近傍配置の場合に生じ、「前後の流れからはじき出される“連行離脱”」
図２の第３行第４列に示す配置の場合に生じるモデルを仮定した。そして、“連行参入”と“連が図行離脱”の確率をいろいろ変化させた。図４から、流体粘性が大きく変化している様子がわかる。

流体解析における時間発展計算の
AI サロゲートモデルの追求― 25 ―

ケース 1
連行参入確率：20.0%
連行離脱確率：10.0％

ケース 4
連行参入確率：28.3%
連行離脱確率：1.7％

ケース 2
連行参入確率：25.0%
連行離脱確率：5.0％

ケース 5
連行参入確率：28.4%
連行離脱確率：1.6％

ケース 3
連行参入確率：28.0%
連行離脱確率：2.0％

ケース 6
連行参入確率：28.6%
連行離脱確率：1.4%

図 ４ ． 負の粘 性 発 現と正 の 粘 性発現 の つ り合い に よ る流体 粘 性 制御の 試 み（衝突散乱後の仮想粒子の“連行参入”が負の粘性を、“連行離脱”が正の粘性を発現させる。）

― 26 ―

SENAC Vol. 55, No. （1 2022. 1）

６．“仮想粒子反転連行モデル”がもつ潜在的なメリットとまとめ（１）理論的な観点
乱流の模擬を実現するにあたって、ナビエ・ストークス方程式を解く通常の数値流体力学(CFD)
による解法との関係が気になるところである。
まず言えることは、Teixeira の“４次元 FCHC54 速度モデル”による格子法流体解析の結果は、
非圧縮性流体に関するナビエ・ストークス方程式に基づく解析結果と比較して、マッハ数に関す[6] 。
る３次の精度で一致することが、理論的に示されている[
さらに、以下のとおり、CFD で必要になるような“複雑な乱流モデル”を考える必要はないと思われる。
“サブマクロ粒子”は、多数のミクロ粒子の集団を仮想的にひとつの巨大な粒子のイメージにまとめたものである。従って、実際には、あるサブマクロ粒子に属するミクロ粒子は、近傍のサブマクロ粒子に属するミクロ粒子と相互作用を起こすので、一般に、サブマクロ粒子がもつ運動量などの物理量は揺らいでいるはずである。原子や分子を想定した粒子に関する物理法則は、通常、このような内部構造に起因した物理量の揺らぎを仮定しないミクロ粒子を適用対象にして設定されたものと考えられる。従って、サブマクロ粒子の物理量についても同じ法則をそのまま適用していいとは限らない。より具体的には、以下のとおりである。
まず、原子や分子のミクロスケールで、物理学的な理論モデルを仮定する。ここでは、格子ガス法を考察しているので、この理論モデルは“ボルツマン方程式”であろう。計算機シミュレー
ションでは、離散化した格子上で計算を行わざるをえないので、実際には、“格子ボルツマン方
程式”をミクロスケールの超微細格子上で計算することが、物理学的理論モデルを忠実に解く方法だと考えられる。従って、本来は、ミクロスケールの超微細格子で“格子ボルツマン方程式”に基づく粒子挙動を計算し、これを合計(平均)してマクロスケールの流体挙動を得る必要がある。
しかし、このような計算は、常識的な時間では実行できない。そこで、格子ガス法では、サブマ
クロ粒子が“格子ボルツマン方程式”を満足すると仮定して、サブマクロスケールから計算を開始する。この結果が、ミクロスケールから時間をかけて計算したと仮定した場合の結果に、理論的に一致するか否か？が問題である。
Teixeira の“４次元 FCHC54 速度モデル”については、サブマクロスケールで平均された仮想粒子の数密度が、ミクロ粒子の数密度が満足すべきと仮定した“格子ボルツマン方程式”を精度よく満足することが示されている[6]。また、このことは、「サブマクロスケールにおいて個々の格子点の物理量が激しく揺らぎ大きなノイズが入っているように見えても、その数値をそのまま合計(平均)して疎視化すれば、必ず適切な結果を得ることができる」ことを意味している。
（２）AI 的な観点
本稿で掲げた具体的な主題は、「NS 方程式による時間発展計算を超高速に実行できる“AI サロ
ゲートモデル”を構築する」ことであった。このため、格子ガス法における仮想粒子のイメージ
から出発することにした。特に、格子ガス法における仮想粒子の衝突散乱と並進移動のイメージ
は、「“仮想粒子が移動するパスをシナプス結合にもつニューロン”が多数結合したニューラルネ
ットワーク」を自然に連想させる。具体的には、格子ガス法の時間発展過程は、リカレント型多[3]。このため、“誤差逆伝搬法”による学習がリアル
層パーセプトロンの時間展開計算に等しい[
タイムデータ同化を実現する基本的手法として利用できる。また、このニューラルネットワークの計算は、入出力と重みの値を微小整数の範囲に限定できるので、演算の簡素化と記憶容量の節約により、高速で低消費電力な計算の実現に大きく寄与できるであろう。さらに、１格子点に関する時間発展計算を１ビット幅で実行できるので超並列計算が可能であることも、従来から述べている[[2-4]。また、ニューラルネットワークの基本構造を、格子ガス法の物理的イメージから導流体解析における時間発展計算の
AI サロゲートモデルの追求― 27 ―

いているので、計測情報を学習して得た新しい確率頻度の物理的な意味を考察できる。すなわち、“説明可能な AI”になるので、次の改良に我々の思考を反映させやすくなる。
（３）将来展望のまとめ
今後の工学的なものづくりをめぐる計算機シミュレーションの将来を展望したとき、「マルチス
ケール、マルチフィジクス、リアルワールド(不完全な計測情報でも機能するリアルタイムデータ
同化等)、カーボンニュートラル(温室効果ガス排出削減のための超低消費電力等)の同時実現」が究極的な達成目標になってくると考えられる。これを実現する道のりは長いかもしれないが、本稿で述べたような“サブマクロ粒子イメージに基づく単純化された AI サロゲートモデル”は、その前進に大きく貢献できると思われる。特に、確率的なビット演算による仮想粒子の衝突散乱計
算は、“ゲート型量子計算機”の応用に適していると思われる。上記の究極的な目標は、結局、量子計算機によって達成されるのかもしれない。
しかし、まずは、足元の第一目標として、幅広いレイノルズ数の流体シミュレーションを可能にするため、“ゼロに近い正粘性を安定的に発現できる方法”が必要であり、この方法を可能な限り単純なモデルで確立し、しっかり検証しておくことが望まれる。ここでは、まだ不十分であるが、その方向のヒントになる“仮想粒子反転連行モデル”を示した。

謝辞
本稿で述べた研究課題の実行にあたっては、“AI サロゲートモデル”の候補となる“モデル”を多数考案し、それぞれのモデルについて、多数のパラメータスタディを行う必要があった。このとき、異なる計算モデルを試すたびに計算コードのチューニングに多くの時間をさくことは非効率的で悩ましい。この点、東北大学サイバーサイエンスセンターのベクトル型スーパーコンピュ
ータ AOBA-A は、特別なチューニングをしないでも十分な計算速度を確保することができ、大変有益であった。また、利用にあたって同センター関係各位のご親切なご指導とご協力をいただき、心から感謝する次第である。今後とも、このような使い勝手のよいベクトル型スーパーコンピュ
ータのさらなる開発導入と同センターの有意義な活動の継続を期待している。
参考文献[1] 松岡,“格子ガス法流体解析モデルとニューラルネットワークの融合”, SENAC Vol.54
No.1, pp.39-49, 2021
[2] 松岡, 菊池,“多速さ格子ガス法実用化展開への手がかり”, SENAC Vol.49 No.4, pp.1-15,
2016[3] 松岡,菊池,“コンパクトな計算機によるリアルタイム流体解析の実現に向けて”, SENAC
Vol.51 No.2, pp.1-10, 2018[4] 松岡,菊池,“仮想粒子の並進移動過程に干渉効果を加味した流体解析の可能性”, SENAC
Vol.52 No.2, pp.18-27, 2019[5] 松岡,菊池,“リカレントニューラルネットワークによる実世界流れ場解析用時間発展計算モデルの探求”, SENAC Vol.53 No.1, pp.25-33, 2020
[6] Chistopher M. Teixeira, “Continuum Limit of Lattice Gas Fluid Dynamics”, MIT,1993
[7] B. Haaslacher, U. Frisch, Y. Pomeau, “Lattice Gas Automata for the Navier-Stokes
Equation”, Physical Review Letters Vol.56, No.14, pp.1505-1508, 1986
[8] Uriel Frisch, Dominique d'Humiè res, Brosl Hasslacher, Pierre Lallemand, Yves
Pomeau,Jean-Pierre Rivet, “Lattice Gas Hydrodynamics in Two and Three Dimensions”,
Complex Systems,１(1987), pp.649-707, 1987
[9] Michel Hénon, “Viscosity of a Lattice Gas”, Complex Systems, 1(1987), pp.763-789,

― 28 ―

SENAC Vol. 55, No. （1 2022. 1）

1987[10] Hudong Chen, Chris Teixeira, Kim Molvig, “Digital physics approach to computational
fluid dynamics: Some basic theoretical features”, International Journal of Modern
Physics C, Vol.8, No.4 (1997), pp.675-684, 1997[11] 松岡,“ビット演算による CFD と等価な高精度流体解析手法”, RIST News No.64, pp.1728, 2018[12]日野,“乱流の科学”, pp597-637, 2020, 朝倉書店（ISBN 978-4-254-20161-1）[13] Daniel H. Rothman, “Negative-Viscosity Lattice Gases”, Journal of Statistical
Physics, Vol.56, Nos. 3/4, 1989
円 柱 後 流 の過渡 変 化 計算に よ る 流体粘 性 の 変化の 観 察 に利用 し た 計算体 系 等【補足】円
粘性制御の効果を確認するため、円柱後流のシミュレーション計算を行った。粘性に変化があれば、後流の様子が変化するはずである。
具体的には、東北大学サイバーサイエンスセンターのベクトル型スーパーコンピュータ AOBA-Aの 16CPU(128 コア)を利用し、MPI による並列計算で約 12 分の小規模計算で求められる円柱後流の３次元過渡変化計算（約 3900 万格子点で 25600 時刻ステップ、円柱軸に垂直なある平面上の流体挙動）を試みた。
計算条件等の詳細を下下図に示す。
３次元空間中で多数の格子点をＸ
ＹＺの各方向に並べ、直方体形状の格Z

子点配列を作る。各格子点は、その内

流体入口

流体出口
X

流れ

部に４次元目の座標として R＝R

Z

0,1,2,3 の位置を識別できる自由度をXもつとする。これが左図であり、４次
元面心超立方体格子を３次元空間へ

Y

Y

Z

投影した姿である。今回の数値シミュ
レーションでは、３次元縮退格子として、Ｘ方向に 1920 個、Ｙ方向に 1280
X個、Ｚ方向に 16 個の格子点を配置した。
格子点間隔の長さを 1 にとると、円柱の直径は 40 に設定した。

Y[ 過 渡変 化シ ミュ レ ーシ ョンの 条 件]
シミュレーション計算を開始する時刻ステップ０の時点で、各格子点には、そこに存在できる仮想粒子の最大数の 20％の数の仮想粒子をランダムな向きの速度で配置する。この結果、疎視化して得られるマクロな流速はゼロであり、流体は、直方体形状の中で静止している。
次に、時刻ステップ１の時点から、＋Ｘ向きの速度をもつ仮想粒子をＸ＝０の位置から注入していく。
すると、時刻ステップが進むにつれて、流体全体が＋Ｘ向きのマクロな速度をもつようになる。このとき、＋Ｘ側の先にある直方体出口においては、出口直前に存在する格子点上の仮想粒子配置を、出口直後に存在する格子点の仮想粒子配置にコピーして、出口におけるマクロな流速の勾配がゼロになるという境界条
件を近似的に実現した。また、±Ｙ方向と±Ｚ方向には、周期的境界条件を適用した。この流れの中の入り口に近い位置に、“Ｚ方向の中心軸をもつ無限大の長さの円柱”を置き、その後流に生じる流体挙動を計算した。

図 ． 円 柱後流 の 過 渡変化 シ ミ ュレー シ ョ ンを行 う 計 算SENAC Vol. 55, No. 1（2022. 1）［大学 ICT 推進協議会― 29 ―

2021 年度年次大会論文集より］

スーパーコンピュータ AOBA の運用状況と利用者支援について
森谷 友映 1) , 佐々木 大輔 1) , 齋藤 敦子 1) , 山下 毅 1) ,小野 敏 1) , 大泉 健治 1) , 滝沢 寛之 2)
1) 東北大学 情報部情報基盤課

2) 東北大学サイバーサイエンスセンター
t-moriya@cc.tohoku.ac.jp

Operation status and user support of Supercomputer “AOBA”
Tomoaki Moriya1) , Daisuke Sasaki1) , Atsuko Saito 1) , Takeshi Yamashita 1) ,Satoshi Ono1) ,Kenji Oizumi1) ,Hiroyuki Takizawa2)
1) Information Infrastructure Division, Information Department, Tohoku Univ.
2) Cyberscience Center, Tohoku Univ.
概要

東北大学サイバーサイエンスセンターは，全国共同利用設備として大規模科学計算システムの整備と，HPCI の資源提供機関としての役割を担っている。本稿では，2020 年 10 月から運用を開始したスーパーコンピュータ AOBAの運用状況と利用者支援について紹介する。

1 スーパーコンピュータ AOBA

東北大学サイバーサイエンスセンター (以下，本置を搭載している。

サブシステム AOBA-A(以下，AOBA-A) は，1VHと 8VE が構成単位となる B401-8 モデルを採用し，センター) では，2020 年 10 月からスーパーコンピ

全体では 72 個の VH と 576 個の VE で構成される。

ュータ AOBA はサブシステム AOBA-A(SX-Aurora

は，1.48PFLOPS(倍精度)，総主記憶容量は 45TB，総ュータ AOBA の運用を開始した [1]。スーパーコンピ

TSUBASA B401-8，日本電気株式会社製)，サブシステ

ム AOBA-B(LX 406Rz-2，日本電気株式会社製) の 2

種類の計算機システムと，ストレージシステム (DDN

SFA7990XE，DDN 社製)，大判プリンタ，講習会端末およびそれらを接続するネットワーク機器群で構成される。図 1 にシステム構成図を示す。

以下では，スーパーコンピュータ AOBA の性能について紹介する。

1.1

サブシステム AOBA-A

VE と VH を合わせたシステム全体の理論演算性能
メモリバンド幅は 895.68TB/s となり，高い演算性能とメモリ性能の両立を実現している。

1.2

サブシステム AOBA-B

LX 406Rz-2 は，1 ノードに AMD EPYC プロセッ

サ 7702(64 コア) を 2 基と 256GB の主記憶装置を搭載し，合計 68 ノードで構成される。ノードあたりの理論演算性能は 4.096TFLOPS(倍精度) ，全体の理論

演算性能は，278.5TFLOPS(倍精度)，総主記憶容量は17TB，総メモリバンド幅は 27.2TB/s となる。サブ

SX-Aurora TSUBASA はアプリケーション演算処システム AOBA-B (以下，AOBA-B) は，ベクトル演処理を行うベクトルホスト (以下，VH) 部により構成

オープンソースソフトウェアの高速な実行を目的とし理を行うベクトルエンジン (以下，VE) 部と，主に OSされる。PCIe カードに搭載される VE 部はベクトルプロセッサおよび高速メモリから構成され，x86/Linuxが動作する VH と PCIe 経由で接続される。

本センターが導入した VE(Type 20B) は，理論演

算性能 2,456GFLOPS(倍精度) となるマルチコア (8コア) ベクトルプロセッサを 1 基と 48GB の主記憶装算に不向きなプログラムや，商用アプリケーション，て導入されている。

2 計算機利用状況
2.1

稼働利用率

スーパーコンピュータ AOBA を導入した 2020 年

10 月から 2021 年 8 月までの月ごとの稼働利用率をサ

― 30 ―

SENAC Vol. 55, No. （1 2022. 1）

図1

スーパーコンピュータ AOBA 構成図

図 2 AOBA-A：稼働利用率の推移

図 3 AOBA-B：稼働利用率の推移

実行時間

ブシステムごとに図 2 および図 3 に示す。ここで稼働

2.2間と，実際に利用された総 VE/ノード時間との比率でジョブの利用形態として以下の 2 つがある。

利用率とは，稼働していた計算資源の総 VE/ノード時あり，サービス提供可能だった計算資源が実際に使われた割合を示している。なお，2020 年 10 月は試験運用期間で算出に必要なデータを採取していなかったた
め，稼働利用率の掲載は割愛する。

図 2 および図 3 から分かる通り，AOBA-A，AOBA-

B 共に 80% を超える高い稼働利用率で運用が続いている。

本センターでは，スーパーコンピュータ AOBA の• 共有利用：利用する VE/ノード数を指定してジョブ投入した際，利用者間で VE/ノードを共有しジョブ実行する。

• 占有利用：VE/ノードを占有して確保し，ジョブ投入をすることができる。その際，他の利用者のジョブ終了を待つ必要がない。

スーパーコンピュータ AOBA の運用状況と利用者支援について― 31 ―

図 6 AOBA-A 実行件数の割合

図 4 AOBA-A 実行件数

図 7 AOBA-B 実行件数の割合

3 利用者支援活動

本センターでは，計算機利用に関する利用相談窓口

図 5 AOBA-B 実行件数

共有利用を対象に，2020 年 10 月から 2021 年 8 月を設けている。相談者は学内外問わず，本センター利用者と利用予定者である。主に，Web フォームまたは事前予約制による面談 (コロナ禍のためオンライン

ツール活用) にて問い合わせを受け付けている。対応までの期間における，実行時間ごとのジョブ件数をサ

は，技術職員とテクニカルアシスタントが行うが，内AOBA-B 共に，利用 VE/ノード数によらず，「24 時

図 8 に，2020 年 10 月から 2021 年 8 月までの問い

ブシステムごとに図 4 および図 5 に示す。AOBA-A，

容によっては，ベンダーへ協力依頼をする。

間以内」で実行終了したジョブが最多であった。

合わせ件数を示す。

れ図 6 および図 7 に示す。「24 時間以内」で実行終了

多く，76 件の問い合わせがあった。多くあった問い合

次に，図 4 および図 5 を比率で表したものをそれぞ

運用開始最初の 2020 年 10 月は集計期間の中で最もしたジョブは実行件数が多く，全体において AOBA-

わせは，コンパイルコマンドの変更，ライブラリの指ジョブの割合を抽出すると，AOBA-B で特徴的な点

ステムからスーパーコンピュータ AOBA へアーキテ

A：99%，AOBA-B：98% を占めた。24 時間を超えるとして，「168 時間以上」のジョブが 12% を占めていた。このことから，シミュレーションを大規模化し，長時間実行する利用者が一定数いることが分かる。

スーパーコンピュータ AOBA では，占有利用，共定方法，ログイン操作に関するものだった。従前のシクチャの変更や，セキュリティ対策の一環としてログ
インノードを設けるなどシステムの構成および利用法に変更があったためと考えられる。

2020 年 11 月と 12 月には，利用者向け講習会を計

有利用どちらでも 1 つのジョブを最大 720 時間まで

7 回開催した。その効果もあり，2021 年 1 月以降，問ジョブの最大実行時間は，AOBA-A では共有利用がトピックは，マニュアルや FAQ 等へフィードバック

実行することができる。そのうち，実際に実行された

360 時間，占有利用が 720 時間，AOBA-B では共有利

用が 440 時間，占有利用が 720 時間であった。

い合わせは落ち着いてきた。今後も問い合わせが多い

していく予定である。

また高速化支援活動として，以下の 2 点を実施して

― 32 ―

SENAC Vol. 55, No. （1 2022. 1）

表1

高速化支援活動実績

直近では 2021 年 3 月と 8 月にそれぞれ実施し，今後も計画予定である。詳細は，8 月時に広報した本セン

ターの大規模科学計算システムの Web ページ*1 を参照いただきたい。

4 まとめ

本稿では，東北大学サイバーサイエンスセンターにおける 2020 年 10 月から 2021 年 8 月までのスーパーコンピュータ AOBA の利用状況を報告し，利用者支援の取り組みについて紹介した。

AOBA-A・AOBA-B 共に稼働利用率では，2020 年

11 月以降 80% を超える利用があり，引き続き維持できるよう努めたいと考える。利用者支援活動では，問図8

問い合わせ件数の推移いる。

1 点目は，利用者のジョブ実行履歴情報から実行性能を抽出し，該当利用者へプログラムの高速化に対するサポートを実施している。これまで本センターが行ってきた，高速化支援の成果を表 1 に示す。1997 年

から 2020 年にかけて 237 件の高速化支援を行った。
単体性能では平均約 15.0 倍の性能向上を，並列性能

い合わせ内容をもとにマニュアルや FAQ のコンテンツを充実させていく予定である。そして，性能向上に関するサポート面では，高速化支援活動を継続して実施していきたい。

参考文献[1] 山下毅，森谷友映，佐々木大輔，齋藤敦子，小野敏，大泉健治，滝沢寛之，「スーパーコンピュータ

AOBA の紹介」，SENAC Vol.54 No.1（2021-1），
pp.50-55，2021．

では約 10.1 倍の性能向上を得ることができ，大幅な向上を実現している。この活動により利用者のプログ

ラムは性能改善しているが，大規模化・長時間化をする様々なプログラムに今後も対応すべく，高速化支援
活動を推進していきたい。

2 点目は，各サブシステムの最大規模（AOBA-A：

576VE，AOBA-B：68 ノード) まで利用できる機会

として大規模並列ジョブチャレンジを実施している。

*1 https://www.ss.cc.tohoku.ac.jp/n20210714-1/

SENAC Vol. 55, No. 1（2022. 1）［大学 ICT 推進協議会― 33 ―

2021 年度年次大会論文集より］

キャンパス外 eduroam と大学における OpenRoaming 導入
後藤英昭 1) ,原田寛之 2) ,漆谷重雄 3)

1) 東北大学 サイバーサイエンスセンター
2) 札幌学院大学 情報処理課
3) 国立情報学研究所

Oﬀ-campus eduroam and Deployment of OpenRoaming atUniversities
Hideaki Goto1) , Hiroyuki Harada2) , Shigeo Urushidani3)
1) Cyberscience Center, Tohoku University
2) Information Processing Division, Sapporo Gakuin University3) National Institute of Informatics
概要

学生や生徒，教職員の学習・教育・研究環境を改善するために，キャンパスのみならず市街地でも eduroam のサービスを提供したいという要求があり，各国が様々な枠組みで対応してきた．これと反対に，他大学からの訪問者に限らず，より広く，市民の学習環境の改善のために，キャンパスに公衆無線 LAN を導入する事例も増えてきている．筆者らは，社会全体の ICT 対応推進という側面を踏まえて，従来型のキャンパス無線 LAN と公衆無線

LAN の間の垣根を減らした，新しい無線 LAN インフラのあり方を提唱してきた．Wireless Broadband Alliance(WBA) により運用が始まった OpenRoaming は，同組織の次世代ホットスポットのコンセプトを具現化した無線
LAN ローミング基盤であり，筆者らもその初期の開発段階から参画している．eduroam と OpenRoaming の連携のために，eduroam の運用母体である GÉANT において相互運用の検討が進められ，試験運用ができる状態になった．本報告では，OpenRoaming 基地局で eduroam を有効にする方法と，キャンパスに OpenRoaming を導入する方法について，最新動向と実践的な技術情報の提供を行う．

1 はじめに教育・研究機関向けの無線 LAN ローミング基盤である eduroam (エデュローム) は，執筆時点 (2021 年

報のとおり，筆者らは，社会全体の ICT 対応推進という側面を踏まえて，従来型のキャンパス無線 LANと公衆無線 LAN の間の垣根を減らした，新しい無線 LAN インフラのあり方を提唱している．Wireless9 月) で世界 106 か国 (地域)，国内 313 機関に導入

Broadband Alliance (WBA) により運用が始まったても，海外では既にいくつかの国で導入が進み，社トのコンセプトを具現化した無線 LAN ローミング基されるに至っている [1]．初等・中等教育機関につい

OpenRoaming [5] は，同組織の次世代ホットスポッ

会全体の学習・教育環境の ICT 対応が進められてい盤であり，筆者らもその初期の開発段階から参画してを改善するために，キャンパスのみならず市街地でも包含する形の，Cityroam [6] と呼ぶローミング基盤を各国が様々な枠組みで対応してきた [3, 4]．この種のる GÉANT において，OpenRoaming との相互運用のと呼ばれるようになった．これと反対に，他大学か報告では，各機関が市街地等の OpenRoaming 基地局

る [2]．学生や生徒，教職員の学習・教育・研究環境

eduroam のサービスを提供したいという要求があり，サービスは，ここ数年でよく “oﬀ-campus eduroam”

いる．また，国内では，eduroam と OpenRoaming を開発・運用している (図 1)．eduroam の運用母体であ検討が進められ，試験運用ができる状態になった．本らの訪問者に限らず，より広く，市民の学習環境の改からの eduroam 接続を受け入れる方法と，自機関の事例も増えてきている．2020 年度年次大会 [2] で既ついて，最新動向と実践的な技術情報の提供を行う．

善のために，キャンパスに公衆無線 LAN を導入するキャンパスに OpenRoaming を導入する方法の二点に

― 34 ―

SENAC Vol. 55, No. （1 2022. 1）

eduroam via OpenRoaming ANPs

Universities providing OpenRoaming
service for visitors and residents.)

IdPs
eduroam ROs
eduroam.org.usDefault
eduroam-OR
gateway

.us

.jp

.jp
via local hubs
or vendors(not necessarily
eduroam’s)

RadSec via hubs
or direct
IdPs

図 2 eduroam と OpenRoaming の認証連携．

のブローカーとなる組織から発行される．

2.2

eduroam の基盤を用いた接続

eduroam の国際的な運用母体である GÉANT と，

国内のセキュア公衆無線 LAN ローミング基盤である
図 1 Cityroam ウェブサイト (2021 年 9 月)．

2 OpenRoaming 基地局からの eduroam 利
用の有効化

2.1

概要

Cityroam は，2020 年春に，OpenRoaming に初期サポートメンバーとして参加した．現在，GÉANT はWBA の Contributor Member であり，OpenRoam-ing においてブローカーの権限も有している．このため，国際的な eduroam としては，そのメンバーに対して RadSec 用電子証明書の発行も可能である．

eduroam メンバーの具体的な OpenRoaming 参加方

本章では，利用者が eduroam のアカウントを用いて法はまだ Global eduroam Governance Committeeの，各機関で必要な対応方法を示す．従来の eduroam

すれば，国単位で大学等の OpenRoaming 利用が可能

IDentiﬁer) が用いられており，この SSID を吹いてい

国内の eduroam 参加機関が，OpenRoaming からOpenRoaming の基地局に接続できるようにするため
では，“eduroam” という共通の SSID (Service Set

る基地局ではすべての eduroam アカウントが利用できるように，eduroam Compliance Statement に規定がある．これに対して，OpenRoaming 上で eduroam

アカウントを使う場合は，ID プロバイダ (IdP) となる各機関で，オプトインが必要である．

eduroam 参 加 機 関 が ，OpenRoaming 上 で のeduroam 利用を有効化するには，OpenRoaming からの認証要求が自機関の ID プロバイダ (IdP) まで届くように，認証経路を設定する必要がある．図 2(GeGC) で検討中であるが，各国のオペレータが希望
な段階にある．

の接続を受け入れたい場合，RadSec エンドポイント

として，以下のような選択肢が考えられる．

• 自機関の proxy を使う．(証明書の取得が必要)• eduroam (国際) の proxy を使う．
• 国内の eduroam proxy を使う．

• その他の事業者の proxy を使う．

現在，国内の eduroam 運用組織である eduroam JP(国立情報学研究所) では，OpenRoaming の扱いが決に，eduroam と OpenRoaming の認証連携の様子をまっていない．また，各機関が証明書を取得するため付ける様子であり，上部の雲の中にある IdP が各機関

が，これには各国の代表組織による取次が必要となる．

に RadSec [7] を用いた Peer-to-Peer 的なルーティン

eduroam JP 宛ての認証要求を中継できる．従って，示す．左図が，OpenRoaming からの認証要求を受けに相当する．OpenRoaming の認証連携では，全面的グが行われるが，各機関の IdP 自身，または，その途中にあるゲートウェイ (プロキシ) が RadSec エンドポ

イントになる必要がある．RadSec エンドポイントにはサーバ証明書が必要であり，これは OpenRoamingには，eduroam の国際運用組織への申請が必要となる

その他の事業者として，現在，国内では Cityroam が現時点で利用できるのは，eduroam (国際) 経由，または，Cityroam 経由の接続のみである．

OpenRoaming 上 で eduroam ア カ ウ ン ト を 有 効に し た い 機 関 は ，RadSec で 利 用 さ れ る NAPTR

キャンパス外 eduroam と大学における OpenRoaming 導入― 35 ―(Naming Authority Pointer) レ コ ー ド を 自 機 関 の

DNS (Domain Name System) サーバに追加する必要がある．OpenRoaming の ANP (Access Network

Provider. eduroam の SP に相当) は，認証要求に含まれるレルムを元に DNS で名前解決を行い，実際にサービスを提供しているサーバのアドレス (ここではプロキシの IP アドレス) を求める．

2.3 Cityroam 経由の接続

執筆時点で，Cityroam は OpenRoaming に接続済

みであり，eduroam JP の参加機関もトライアルが可能である．Cityroam を中継点とする場合，具体的には，以下のようなレコードを登録すればよい．(機関側図3

新宿中央公園の Cityroam 対応スマートポール (東京都デジタルサービス局)．

で SRV レコードの登録は不要)<レルム> IN NAPTR 50 50 "s""aaa+auth:radius.tls.tcp" ""_radiustls._tcp.idp.cityroam.jp.

どのようにしてプロキシのアドレスが取得されるのか，dig コマンドを用いて具体例を見ることができる．

$ dig naptr cityroam.jp$ dig srv _radiustls._tcp.jpgw.cityroam.jp
Cityroam 経由の接続を行う場合，事前の報告は不要であるが，長期間利用にはご一報いただきたい．

2.4

接続用プロファイル

OpenRoaming の利用には，Passpoint [8] に対応し

図4

札幌学院大学における OpenRoaming．

れる．

図 3 は，東京都が実証実験中のスマートポールの例

で，eduroam と Cityroam に加えて，OpenRoaming

た端末と，Passpoint プロファイルが必要である．各機にも対応した基地局が内蔵されている．国際ローミンのは難しいと考えられるが，幸いなことに，eduroam

国内外で普及することを願っている．

関でこのプロファイルを発行するシステムを用意する

CAT (Conﬁguration Assistant Tool) が既に対応しており，従来の eduroam のプロファイルに加えて，OpenRoaming 用のプロファイルも端末に同時インス

トール可能となっている．

2.5

サービスエリア

グに対応したこのような公衆無線 LAN サービスが，3 自機関への OpenRoaming 導入

大学等が OpenRoaming 対応の基地局を導入して市民向けの公衆無線 LAN サービスを提供する方法を紹介する．図 2 の右側に，機関が OpenRoaming

Cityroam では eduroam の SSID も併設することをANP となる場合の認証連携の様子を示す．左図のing 経由の接続を試すには，eduroam SSID 用のプロの雲は eduroam ではなく，OpenRoaming のもので用については eduroam SSID の方が簡便なため，わOpenRoaming に参加する必要がある．

基本としているため，Cityroam 基地局で OpenRoam-

ファイルを削除する必要がある．もちろん，実際の利

ざわざ OpenRoaming 経由で利用する必要はないと言える．

国際的にも OpenRoaming は黎明期にあり，サービ

上部の雲は eduroam の基盤を表しているが，右上ある．OpenRoaming のサービスを提供する機関は，前章で述べたように，eduroam は全体として既にOpenRoaming に参加している．しかしながら，各国のオペレータは eduroam の運用組織であり，市民一スエリアが少ない．今後，観光の復活に伴って徐々に般を受け入れるような OpenRoaming 基地局の運用はも，Cityroam 以外の基地局が増えてくれば，Open-OpenRoaming の認証要求をそのまま重畳させること

対応基地局が増えることが期待されている．国内でRoaming 経由の接続が有用になってくると予想さ管轄外である．eduroam の認証連携ネットワークに，は，技術的に困難である．

― 36 ―

現在，大学等が OpenRoaming の基地局を導入するには，以下のような方法がある．

• ハブ事業者経由で OpenRoaming の認証連携基盤に接続する．

• OpenRoaming に対応した無線 LAN 事業者 (通信事業者) に基地局システムの運用を委託する．

• 機関自身が OpenRoaming に参加する．

• Cisco DNA Spaces のような，OpenRoaming が統合された無線 LAN システムを導入する．

機関自身がブローカーを通じて OpenRoaming に参加する方法がある．しかし，現時点でこれに対応できる

国内のブローカーは知られていない．ブローカーとしての eduroam から RadSec の証明書を取得することも考えられるが，現時点ではまだ eduroam に運用の枠組みがなく，計画も未定である．

も し ，機 関 自 身 が 運 用 し て い る 基 地 局 に Open-

Roaming のサービスを追加するのであれば，City-

roam 経由で OpenRoaming に接続するのが最も簡便と考えられる．この場合，Cityroam がハブ事業者として機能する．自前で基地局を運用する必要がないのであれば，通信事業者に業務委託する方法が，安定性でも管理・運用面でもメリットがあると考えられる．

図 4 は ，札 幌 学 院 大 学 の キ ャ ン パ ス に 導 入 し た

OpenRoaming の例である．認証連携には，Cityroamの基盤を利用している．OpenRoaming では様々な機

関のアカウントが利用できるが，最近よく利用されているのは，Google や Apple，Samsung のアカウントを用いてサインアップできる，OpenRoaming 用アカ
ウントである．

OpenRoaming に使うバックホール回線には，市民のトラフィックを流せるものが必要となる．詳しくは2020 年度の報告 [2] をご覧いただきたい．

4 むすび

本報告では，OpenRoaming 基地局で eduroam を有効にする方法と，キャンパスに OpenRoaming を導入する方法について，最新動向と実践的な技術情報を示した．市民向けの国際無線 LAN ローミング基

SENAC Vol. 55, No. （1 2022. 1）

必要がある．筆者らが開発・運用している Cityroam

基盤では，eduroam と OpenRoaming を統合した認証連携基盤を提供しており，eduroam JP の参加機関

であれば，OpenRoaming の試行も可能である．

地域社会への貢献や災害対応などの観点で，学校や図書館等公共施設への安全な公衆無線 LAN の導入が望まれている．今後，国内外の通信事業者とのローミングを充実させ，利用できるアカウントとサービスエ
リアの拡充を図っていく予定である．

本研究の一部は，令和 3 年度国立情報学研究所公募

型共同研究の助成を受けた．

参考文献[1] eduroam JP: https://www.eduroam.jp/(2021 年 9 月 27 日参照)

[2] 後藤英昭, 原田寛之, 漆谷重雄, “キャンパス無線 LAN と公衆無線 LAN の統合 – eduroam とCityroam, OpenRo aming –,” 大学 ICT 推進協
議会 2020 年度年次大会 論文集 WA2-1, 2020.

[3] GÉANT NETWORKS, “eduroam:

Bring-

ing City-Wide Wi-Fi Roaming Access toCape Town.” https://impact.geant.org/
portfolio/eduroam_capetown/(2021 年 9 月 27 日参照)

[4] eduroam Everywhere:

https://www.heanet.ie/services/connectivity/eduroam-everywhere(2021 年 9 月 27 日参照)

[5] WBA OpenRoaming:

https://wballiance.com/openroaming/(2021 年 9 月 27 日参照)

[6] Cityroam: https://cityroam.jp/(2021 年 9 月 27 日参照)

[7] S. Winter and M. McCauley,“DynamicPeer Discovery for RADIUS/TLS and RADIUS/DTLS Based on the Network Access Identiﬁer (NAI),” IETF RFC7585, 2015.
https://tools.ietf.org/html/rfc7585(2021 年 9 月 27 日参照)

盤である OpenRoaming を用いると，eduroam のア[8] Wi-Fi Alliance, “Passpoint – Wi-Fi ホットス

LAN 利用が可能となる．これを実現するためには，な接続を実現.” https://www.wi-fi.org/ja/

カウントを用いて，市街地等の基地局で安全な無線機関ごとのオプトインが必要である．キャンパス内にOpenRoaming 対応の基地局を導入する場合，各機関

は OpenRoaming の認証連携ネットワークに参加する

ポットネットワークへのシームレスでセキュア

discover-wi-fi/passpoint/(2021 年 9 月 27 日参照)

SENAC Vol. 55, No. 1（2022. 1）［報― 37 ―

告］第 32 回高性能シミュレーションに関するワークショップ(WSSP)を開催しました
スーパーコンピューティング研究部 滝沢寛之
東北大学サイバーサイエンスセンターは、ドイツのシュトゥットガルト大学高性能計算センター(HLRS)、学際大規模情報基盤共同利用・共同研究拠点(JHPCN)、HPCI コンソーシアムおよび NEC のご協力を得て、2021 年 12 月 13 日(月)～14 日(火)に高性能計算に関する国際ワークショップ「第 32 回 Workshop on Sustained Simulation Performance (WSSP)」を開催しました。本ワークショップは、国際的に活躍している計算科学の研究者及びスーパーコンピュータ
設計者を招いて、高性能・高効率大規模科学計算に関する最新の研究成果の情報交換を行うとともに、今後のスーパーコンピュータの研究開発のあり方を議論することを目的としています。
今回の WSSP は、東北大学サイバーサイエンスセンターの関係者と日本国内の講演者の一部のみが仙台会場に集まり、それ以外の講演者や聴講者はオンライン参加するというハイブリッ
ド形式での開催となりました。ヨーロッパとの時差を考えて、双方の参加者が比較的参加しやすい時間帯に講演時間を設定しました。日本標準時での開催時間は以下の通りです。
12/12(月) 15：00 - 19：30
12/13(火) 15：00 - 20：30
本ワークショップでは技術講演として全部で 17 件の発表があり、日本、ドイツ、ロシアの研究者により、HPC 技術動向、HPC システム評価、アプリケーション開発の幅広い分野のトピックの講演がありました。参加登録者数は 106 名に上り、終始多くの方々にご参加いただけました。
海外からは、HLRS センター長の Michael Resch 氏による HLRS の現状と課題、および将来展
望に関する講演、モスクワ州立大学研究計算センター長の Vladimir Voevodin 氏による地球物
理学に基づく大規模気象シミュレーションに関する講演、アーヘン工科大学の Matthias
Meinke 氏による飛行機着陸脚の騒音解析に関する講演などがありました。
日本からは、理学研究所の佐藤三久氏、坪倉誠氏、佐野健太郎氏、筑波大学の朴泰佑氏、神戸大学の横川三津夫氏、東京理科大学の塚原隆裕氏、大阪大学の伊達進氏に加えて、多数の NECの技術者が講演を行いました。さらに、主催である東北大学サイバーサイエンスセンターからも、スーパーコンピュータ AOBA の運用や性能評価に関する講演を行いました。
コロナウイルス感染拡大防止の観点から一部の講演者・関係者だけになってしまいましたが、今回の WSSP では仙台会場に集まって議論を深めることができ、実際に会って話をすることの重要性を改めて実感する機会となりました。一刻も早く世界がコロナ前の状況に戻り、以前のような国際会議・ワークショップを開催できるようになることを願うばかりです。
第 32 回 WSSP に関するその他の詳細は、以下のページをご覧ください。
https://www.sc.cc.tohoku.ac.jp/wssp32/ja/index.html

― 38 ―

SENAC Vol. 55, No. （1 2022. 1）

SENAC Vol. 55, No. 1（2022. 1）― 39 ―― SENAC 執筆要項 ―１．お寄せいただきたい投稿内容
サイバーサイエンスセンターでは、研究者・技術者・学生等の方々からの原稿を募集しております。
以下の内容で募集しておりますので、皆さまのご投稿をお待ちしております。なお、一般投稿いただ
いた方には、謝礼として負担金の一部を免除いたします。
・一般利用者の方々が関心をもたれる事項に関する論説・センターの計算機を利用して行った研究論文の概要
・プログラミングの実例と解説・センターに対する意見、要望
・利用者相互の情報交換

２．執筆にあたってご注意いただく事項(1)原稿は横書きです。
(2)術語以外は、「常用漢字」を用い、かなは「現代かなづかい」を用いるものとします。
(3)学術あるいは技術に関する原稿の場合、200 字～400 字程度のアブストラクトをつけてください。
(4)参考文献は通し番号を付し末尾に一括記載し、本文中の該当箇所に引用番号を記入ください。
・雑誌：著者,タイトル,雑誌名,巻,号,ページ,発行年
・書籍：著者,書名,ページ,発行所,発行年
３．原稿の提出方法
原稿のファイル形式は Word を標準としますが、PDF での提出も可能です。サイズ*は以下を参照してください。ファイルは電子メールで提出してください。
－用紙サイズ・文字サイズ等の目安－
・サイズ：A4
・余白：上＝30mm 下＝25mm 左右＝25mm 綴じ代＝0
・標準の文字数（45 文字 47 行）
・表題＝ゴシック体 14pt 中央 ・副題＝明朝体 12pt 中央
・氏名＝明朝体 10.5pt 中央
・所属＝明朝体 10.5pt 中央
・本文＝明朝体 10.5pt
・章・見出し番号＝ゴシック体 11pt～12pt*余白サイズ、文字数、文字サイズは目安とお考えください。
４．その他(1)一般投稿を頂いた方には謝礼として、負担金の一部を免除いたします。免除額は概ね 1 ページ
1 万円を目安とします。詳細は共同利用支援係までお問い合わせください。
(2)投稿予定の原稿が 15 ページを超す場合は共同利用支援係まで前もってご連絡ください。
(3)初回の校正は、執筆者が行って、誤植の防止をはかるものとします。
(4)原稿の提出先は次のとおりです。
東北大学サイバーサイエンスセンター内 情報部情報基盤課共同利用支援係
e-mail uketuke@cc.tohoku.ac.jp
TEL022-795-3406

― 40 ―

SENAC Vol. 55, No. （1 2022. 1）

スタッフ便り2021 年 12 月 1 日付でスーパーコンピューティング研究部に着任いたしました。これからどうぞよろしくお願いいたします。
ここでは、私の趣味の 1 つである料理についてお話ししたいと思います。なぜ料理が好きなのか考えたとき、美味しい食事が好きだから、というのはもちろんなのですが、もう一つの理由としてプログラミングと似ているから、ということがあります。レシピ（アルゴリズム）があり、複数の工程を同時に進めたり（並列処理）、既製品を活用したり（ライブラリ）、プログラミングとの共通点は多いように感じます。事実、表紙に動物の絵が描かれている技術書で有名な出版社は、「プログラマ向けの料理本」を出版しています。
最近は、動画共有サイトでレシピ動画を観て作ることが多くなりました。プロの料理人の動画
は、材料や手順を知るだけでなく、昔なら弟子入りしなければ知りようがなかったような、細かなテクニックを学べるため、勉強になります。プロの調理を見ていると、家庭料理と比べ根本的な材料や工程の違いがあるわけでなく、ちょっとした材料へのこだわりや、調理の手間をかけていることに気づきます。そのような小さな工夫や手間が積み重なり、「プロの味」を作るのではないでしょうか。私は料理人ではないので料理にそこまでの労力を費やせないことも多いですが、業務においては「プロの味」を出せるよう励みたいと思います。 (K.T)
コロナ禍になって 2 回目の新年を迎えました。テレワーク、オンライン会議など仕事の中で新しい経験もありました。個人の生活では遠出の機会が減る中、週末の散歩が恒例となりました。近所に住む友人が 9:50 に迎えに来ます。
なぜ 9:50 かといえば、近所にあるクロワッサン専門店が 10:00に開店するからです。2 人で足取りも軽く 5 分くらいでお店の前に到着。ほとんどの場合先客が店の前にいます。私達の目的は焼き立てのプレーンクロワッサンです。いつも 1 人 1 個しか購入しないのに、お店の方は、毎回爽やかな笑顔と楽しい会話で対応してくれて、それも買い続ける理
由です。県外ナンバーの車や、1 人で 1 万円以上購入する方をみかけることもあります。散歩の前置きが長くなりましたが、それを持って店から 10 分くらい歩いた先の公園を 1,2 周、東屋で持参したコーヒーとともにサックサクのクロワッサンを頂きます。1 年通して歩くと自然の変化がよくわかり、特にカルガモの親
子、季節の花、秋のドングリ拾いは楽しみです。ちっとも痩せない散歩ですが、コロナがなければ知らずに過ごした楽しみを１つみつけました。(K.S)

サイバーサイエンスセンタースタッフ着任のお知らせ
2021.12.1 付け[着任]
髙橋 慧智 スーパーコンピューティング研究部助教（奈良先端科学技術大学院大学助教から）SENAC 編集部会
滝沢寛之 水木敬明 後藤英昭 伊藤昭彦
早坂和勝 大泉健治 小野 敏 斉藤くみ子

編集・発行

PDF 作成令和 4 年 1 月発行
東北大学
サイバーサイエンスセンター
仙台市青葉区荒巻字青葉 6-3
郵便番号 980-8578
株式会社 東誠社

ࢫ࣮ࣃ࣮ࢥࣥࣆ࣮ࣗࢱ $2%$ ࢩࢫࢸ୍࣒ぴィ⟬ᶵࢩࢫࢸ࣒ᶵ ✀ࢧࣈࢩࢫࢸ࣒ $2%$$

6;$XURUD768%$6$

ࢧࣈࢩࢫࢸ࣒ $2%$%/;5]

ࢧ࣮ࣂ࡜࣍ࢫࢺྡ
ࣟࢢ࢖ࣥࢧ࣮ࣂ

ORJLQFFWRKRNXDFMSࢹ࣮ࢱ㌿㏦ࢧ࣮ࣂ

ILOHFFWRKRNXDFMSࢧ࣮ࣅࢫ᫬㛫฼⏝ࢩࢫࢸ࣒ྡ➼฼⏝᫬㛫ᖏ

ࢧࣈࢩࢫࢸ࣒ $2%$$

㐃

⥆

㐠

㌿

ࢧࣈࢩࢫࢸ࣒ $2%$%

㐃

⥆

㐠

㌿

ࢧ ࣮ ࣂ

㐃

⥆

㐠

㌿

㤋ෆ฼⏝

ᖹ᪥ 㹼

ࢧࣈࢩࢫࢸ࣒ $2%$$ ࡢ฼⏝ᙧែ࡜ไ㝈್


฼⏝ᙧែ
↓ᩱ

࣮࢟ࣗྡ 9( ᩘͤ 
V[I



9(



9(

V[
ࠥ

ඹ᭷V[PL[

ᐇ⾜ᙧែ

ࠥ

༨᭷

᭱኱⤒㐣᫬㛫᪤ᐃ್᭱኱್

࣓ࣔࣜࢧ࢖ࢬ ᫬㛫 ᫬㛫9( ༢఩࡛☜ಖ
㸦9+ ࢆඹ⏝ࡋ࡞࠸㸧  ᫬㛫 ᫬㛫 *%9( ᩘ9( ༢఩࡛☜ಖ
㸦9+ ࢆඹ⏝ࡍࡿ㸧ಶูタᐃͤ 9(௨ୖࢆ฼⏝ࡋࡓ୪ิᐇ⾜࡟ࡣ03,ࡢ฼⏝ࡀᚲ⏝ࢧࣈࢩࢫࢸ࣒ $2%$% ࡢ฼⏝ᙧែ࡜ไ㝈್
฼⏝ᙧែ

࣮࢟ࣗྡ

ඹ᭷

O[

ࣀ࣮ࢻᩘͤ ࠥ᭱኱⤒㐣᫬㛫᪤ᐃ್᭱኱್

࣓ࣔࣜࢧ࢖ࢬ ᫬㛫 ᫬㛫*%ࣀ࣮ࢻᩘ

༨᭷

ಶูタᐃ

ͤ ࣀ࣮ࢻ௨ୖࢆ฼⏝ࡋࡓ୪ิᐇ⾜࡟ࡣ03,ࡢ฼⏝ࡀᚲ⏝

Vol.55 No.1 2022-1［巻頭言］
サイバーサイエンスセンターの改組について………………………………菅沼

拓夫

1

強

3

貴志

10

浩

15

敦子

29

［共同研究成果］
ナノアンテナの数値解析のためのブロックモーメント法の高速化…………………………………………………………今野佳祐・陳

細線化ならびに多層化されたプリント基板の FDTD 解析における高速化の検討…………………………………………………………………………春日流体解析における時間発展計算のＡＩサロゲートモデルの追求―直感による究極的な単純化はどこまで有効か？― ………………………松岡［大学 ICT 推進協議会

2021 年度年次大会論文集より］

スーパーコンピュータ AOBA の運用状況と利用者支援について…………………………………………森谷山下
滝沢友映・佐々木大輔・齋藤
毅・小野

敏・大泉

健治

寛之・漆谷

重雄

寛之

キャンパス外 eduroam と大学における OpenRoaming 導入…………………………………………後藤
［報

英昭・原田

33

告］

第 32 回高性能シミュレーションに関するワークショップ（WSSP）を開催しました…………………………………………………………………………滝沢寛之

37

執筆要項 …………………………………………………………………………………………

39

スタッフ便り ……………………………………………………………………………………

40