学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月

jh150030-NA18

階層分割型数値計算フレームワークを用いた 3 次元電磁界解析の高速化研究
杉本振一郎（諏訪東京理科大学工学部機械工学科）
概要 本研究では，演算カーネルに係数が複素数となる疎行列ベクトル積を持つ並列電磁界解析ソルバ ADVENTURE_Magnetic の高速化を行うため，問題規模・並列台数を変えた数値実験を通して大規模電磁界解析の性能向上に向けた問題点の洗い出しを行い，それらを改善することで最終的にピーク性能比 10%程度を目指す．計画に沿って，これまでにウィークスケーリングによる性能評価と複素数演算のチューニングを行った．ウィークスケーリングは FX10 上で最大 1,440 ノードを用いて評価を行い，比較的新しい
複素数演算の実装では 360 ノード以上で性能が低下し，通信性能に課題があることが分かった．今後は非同期通信の導入などにより性能向上を図る．複素数演算のチューニン
グは FX10 と FX100 にて 4 ノード，400 万自由度規模で行い，チューニング前に比べて最大で 1.80 倍の高速化を達成した．しかし FX100 でのピーク性能比は FX10 の半分以下であり，まだ課題が残っている．また別の検討として複素数演算のチューニングを
FX10 8 ノード，100 万自由度規模で行い，ピーク性能比を 7.41%まで向上させることができた．今後はさらに SIMD 演算の幅が倍になった FX100 にも対応するためさらに
SIMD 化を進める検討を行うとともに，より高並列な環境での性能向上を検討し，ピー
ク性能比 10%を目指す．

1.

タメディア工学科）：連続体力学向け DSL の共同研究に関する情報(1) 共同研究を実施した拠点名

開発，高速化東京大学

武居周（宮崎大学工学教育研究部電気システム

名古屋大学

工学科）
：物理モデルの開発，並列化性能評価
室谷浩平（東京大学大学院工学系研究科システ(2) 共同研究分野


ム創成学専攻）：大規模プリプロセス支援

超大規模数値計算系応用分野

片桐孝洋（東京大学情報基盤センタースーパー
コンピューティング研究部門）：高速化(3) 参加研究者の役割分担

2.

杉本振一郎（諏訪東京理科大学工学部機械工学研究の目的と意義

科）：全体統括，並列化コーディング，物理モ(ア) 研究の目的

デル開発

近年，変圧器や回転機などの各種産業用機器，

荻野正雄（名古屋大学情報基盤センター大型計

電車内などの環境中における電磁界解析を効率

算支援環境研究部門）：階層領域分割型ソルバ

的に行うために有限要素法が広く用いられるよーライブラリの開発，並列化性能評価

うになってきた．しかし，電磁界解析分野はそ塩谷隆二（東洋大学総合情報学部総合情報学の問題の複雑さから，構造解析分野などに比べ科）：階層領域分割型ソルバーライブラリの開て分散並列化や高速化の研究が十分ではなく，発，連続体力学向け DSL の開発

インフラとして整備されつつあるスーパーコン

金山寛（日本女子大学理学部数物科学科）：物理

ピュータを活用しきれていないのが現状である．
よって，有限要素法による 3 次元電磁界解析のモデルの開発HPC 利用技術が確立されることは，特に産業界

河合浩志（諏訪東京理科大学工学部コンピュー1

学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月において高く期待されている．そこで我々は大事前の数値実験でその演算性能は計算機環境に

規模数値計算システムの基盤技術として，塩大きく依存することが分かっており，最適化な

谷・荻野らによって有限要素法による構造解析

らびに効率的なコーディング方法の確立が必要分野向けに開発されてきた階層型領域分割法のである．さらに，並列化効率は階層分割型数値技術を応用した，階層分割型数値計算フレーム

計算フレームワークによって十分な性能が得ら

ワークの研究開発を行っている．階層型領域分れることが期待されるが，Tofu などネットワー

割法は Glowinski らが提案した領域分割法を分クトポロジに対する性能評価と最適化は行われ

散メモリ型並列計算機環境に効率よく実装するていない．本研究では，問題規模・並列台数を手法であり，大規模問題を高い並列効率で数値変えた数値実験を通して大規模電磁界解析の性計算できる手法として知られている．本研究で能向上に向けた問題点の洗い出しを行い，それは，主に通常の節点要素や実数向けに開発されらを改善することで，最終的にピーク性能比てきた階層分割型数値計算フレームワークを，10%程度を目指す．

電磁界解析などに出てくる辺要素や複素数を用いた数値計算にも対応させ，3 次元電磁界解析向(イ) 研究の意義

けの並列有限要素解析ソルバとして整備し，ソ

電磁界解析では電場，磁場に空間的な広がり

ルバ全体での演算効率・並列効率を向上させるがあるため，解析したい機器だけでなくその周ことを目標とする．

辺の空間も解析対象としなければならず，解析

我々は階層分割型数値計算フレームワークを領域が大きくなりやすい．また電磁波を対象と用いた並列電磁界解析ソルバ ADVENTURE_

する場合には要素の辺長を電磁波の波長の 10 分

Magnetic (以下，AdvMag)の研究開発を行っての 1 以下にする必要があり，大きく広がる空間

いる．これまでに東京大学情報基盤センター

を大きな要素で埋めることができない．さらにFujitsu PRIMEHPC FX10 (以下，FX10)にて 35

電磁場の様子を詳細に知りたい箇所は機器の中億複素自由度の時間調和渦電流解析，2.6 億複素でも特に薄い構造をしているうえ，電磁界解析

自由度の高周波電磁波解析に成功しているが，で用いる反復法の収束性はアスペクト比に非常

ソルバ全体のピーク性能比は 1～2%程度と計算に敏感であるため扁平な要素でそれらをモデル

資源を有効活用できていない状況である．これ化することができない．そのため電磁界解析のらの解析では本来非定常である問題を，商用電ためのメッシュは大自由度になりやすく，機器

源および電磁波が正弦波的に規則正しく変化すの丸ごと解析では簡単に数千万～数億自由度のることを利用して，複素数での求解を一度だけモデルとなってしまう．しかし，計算時間や大行う準定常問題としている．最終的に得られる規模問題における収束解を得ることの難しさか

連立一次方程式は共役直交共役勾配法などの反ら電磁界解析での並列化研究はなかなか浸透せ復解法で解くことができるため，数値計算におず，企業における設計・開発の現場ではいまでける演算カーネルは係数が複素数となる疎行列

も数千～数十万自由度での解析が主流である．

ベクトル積となる．節点要素を用いた有限要素

そのような自由度に抑えるためにモデルの簡略法における疎行列ベクトル積の高速化に関する化や扱う物理現象の単純化が行われ，十分とは研究は行われているが，辺要素を用いた場合の言いがたい解析精度で工夫を凝らして設計・開性能評価や高速化技術の開発は行われていない．

発が行われているのが現状である．

また，AdvMag は C 言語における複素数演算の本研究によって AdvMag の演算効率，並列効

コーディングとして関数，マクロ，C99 の複素率を向上させることにより，数億自由度の機器

数ライブラリにもとづく実装を持っているが，の丸ごと解析をモデルの簡略化や物理現象の単2

学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月

純化をせずに行うことがより容易になり，大規

節点数，部分領域数を示す．部分領域あたりの模電磁界解析を浸透させる一助になると考える．

要素数が約 100 となるように領域分割を行い，ノードに 2,000 ずつ部分領域を割り当てること

3. 当拠点公募型共同研究として実施した意義

でノードごとの負荷を均一にしている．

物理モデル開発研究が進行中とも言える電磁界解析分野が日々進化する計算機アーキテクチ
ャに対応するためには，工学・応用数学者らと計算科学者らの学際的共同研究体制が必要不可
欠である．本研究体制では，物理モデル開発と並列プログラム実装に強い研究代表者(杉本)を筆頭に，基盤とする階層分割型数値計算フレー
ムワークの開発者であり並列計算に強い研究者

図 1. 無限長ソレノイドコイル(荻野・塩谷)，物理モデル開発の研究者(金山・武居)，高速化の研究者(河合)が参加している．
また，電磁界解析分野はメッシュ品質が非常に重要となるため，プリポストに強い研究者(室谷)も参加している．さらに，複素数演算を高性能かつ効率的に実装するために自動チューニング
技術の研究者(片桐)も参加している．これにより，モデリングと数値計算アルゴリズムの両方を改善しながらシミュレーション全体を高速化でき

図 2. ケーキモデル

た．

表 1. 要素数，節点数，部分領域数

4. 前年度までに得られた研究成果の概要

Comp.
nodes

Elements

Nodes

Subdomains

45

8,739,045

11,733,539

45×2,00090

17,532,512

23,453,016

90×2,000(ア) 性能評価

180

34,665,506

46,232,824

180×2,000

まず FX10 を用いて最大 1,440 ノードのウィ360

69,310,254

92,212,658

360×2,000

720

140,260,096

187,852,707

720×2,000

1,440

277,324,048

370,872,163

1,440×2,000

今年度新規課題のため，該当項目なし．

5. 今年度の研究成果の詳細

ークスケーリングを行った．
モデルとして図 1 の無限長ソレノイドコイルを用いた渦電流解析の精度検証問題[1]を用いる．
導体部の半径は 0.1[m]であるとする．磁気抵抗

インターフェース問題には COCG 法に基づく

は解析領域全体で(1/(4 ))×107 [m/H]，導体

階層型領域分割法を適用し，近似対角スケーリ

は 7.7×106 [S/m]，角周波数はング前処理[2]を用いる．収束判定値は与えず，

2 ×60 [rad/s]とする．コイルに流れる強制電流

反復を 10,000 回繰り返すこととする．部分領域の実部，虚部の大きさは 50，0 [A/m2]と解法には ICCOCG 法を用い，収束判定値を 109，

する．問題の対称性を考慮し，中心角 20°,高さ

前 処 理 の 加 速 係 数 を 1.2 と す る ． 使 用 す る

0.1 [m]の領域を解析対象のモデルとする(図 2

AdvMag のバージョンは 1.6.0 であり，モジュー参照)．表 1 に本報告で用いるメッシュの要素数，

ルはハイブリッド並列版を用いる．また，部分

率

部の導電率
密度

3

学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月

領域間の自由度のばらつきを抑えるため，不導typedef double complex Complex ;体領域にも導電率 0 [S/m]を与える．
#define ComplexPlus( __x, __a, __b ) { ¥struct Complex {*__x = __a + __b ; ¥double re ; /* real part */}double im ; /* imaginary part */};#define ComplexMulti( __x, __a, __b ) { ¥*__x = __a * __b ; ¥図 3. 構造体“Complex”}

図 6. double complex 型とマクロ(C99&Macros)

void ComplexPlus( Complex *x, Complex a,Complex b )C 言語には C99 規格(ISO/IEC 9899:1999)で{(*x).re = a.re + b.re ;標準化された複素数型，例えば倍精度複素数型(*x).im = a.im + b.im ;として double complex 型があり，複素数を実装するにはこの型を用いることが考えられる．し} /* addition */かし，AdvMag の開発が開始されたのは 1990 年
void ComplexMulti( Complex *x, Complex a,

代であり，その当時はまだ C99 規格は存在しな

Complex b )かった．そのため，AdvMag では倍精度の浮動{小数点型 double を 2 つもつ構造体“Complex”(*x).re = a.re*b.re - a.im*b.im ;を用いて複素数を表現した(図 3)．この構造体で(*x).im = a.im*b.re + a.re*b.im ;はメンバ“re”が複素数の実部，“im”が複素数の虚部を表す．また，構造体を用いるために} /* multiplication */AdvMag では四則演算を始め複素数の演算はす

図 4. 構造体と関数(Str.&Func.)

べて関数によって表現していた(Str.&Func.，図4)．しかし関数呼び出しではコンパイラによる最#define ComplexPlus( __x, __a, __b ) { ¥適化が効きにくいため，これらの関数をマクロ(*__x).re = __a.re + __b.re ; ¥によってインライン展開した(Str.&Macros，図(*__x).im = __a.im + __b.im ; ¥5)．さらに構造体を double complex 型に置き換}えた実装も行った(C99&Macros，図 6)．これらの実装と計測は FX10 上で 2014 年度までに行い，#define ComplexMulti( __x, __a, __b ) { ¥ピーク性能比はインライン展開によって 4 倍，

double __re = ¥double complex 型の使用によりさらに 1.15 倍と__a.re*__b.re - __a.im*__b.im ; ¥なった．本報告では Str.&Func.，Str.&Macros，(*__x).im = ¥C99&Macros それぞれについてウィークスケー__a.im*__b.re + __a.re*__b.im ; ¥リングを行う．

(*__x).re = __re ; ¥図 7 にウィークスケーリングの結果を示す．

}

Str.&Func.では 1,440 ノードまで良好なスケー

図 5. 構造体とマクロ(Str.&Macros)

リング結果が得られているが，関数をマクロで
イ ン ラ イ ン 展 開 し た Str.&Macros とC99&Macros は 360 ノード以上で大きく低下し4

学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月

1.2
1

Scalability

0.8
Str.&Func.

0.6

Str.&Macros
C99&Macros

0.4
0.2
0
23

45

90

180

360

720

1,440

2,880

Computation nodes

図 7. ウィークスケーリングの結果d)ている．これは単体性能が低く全体の計算時間c)+スレッド内でのメモリ確保／解放コストに占める通信関連部分の割合が低かったを削減するために実行時に以下を指定

Str.&Func.に対して，インライン展開によって

XOS_MMM_L_ARENA_LOCK_TYPE=0

単体性能が向上したため，通信関連部分が全体

表 2 にチューニング結果を示す．FX10 ではに占める割合が向上したことにより性能低下がc)まででピーク性能比が 1.42 倍，FX100 では d)顕在化したためと思われる．今後は，非同期通までで 1.80 倍になっている．FX100 では FX10

信の導入などにより，性能の改善を図っていく．

に比べてピーク性能比が半分以下であるため，
FX10 から FX100 へのピーク性能の向上に比べ(イ) 複素数演算のチューニング 1て，計算時間の短縮幅は小さいものにとどまっ

FX10 および，名古屋大学情報基盤センターのた．

FX100 にて，複素数演算のチューニングを行った．使用する AdvMag のバージョンは 1.6.0 で

表 2. チューニング結果 1あり，モジュールはハイブリッド並列版を用い

FX10
Time [s] FLOPS /(ratio toPEAK
FX100)[%]
424
4.22(2.27)
412
4.35(2.28)
320
5.57(2.50)
298
5.98(2.37)-

る．複素数の実装は C99&Macros とし，コンパ
イラオプションとして-Kfast を用いてコンパイ

asis

ルしたものを基本とする．モデルとして 400 万
自由度のケーキモデルを用い，部分領域数を

a)30,720 とする．計測には FX10，FX100 とも 4b)

ノードを用いる．

c)

本項での検討として，以下を行う．
a)
b)
c)

d)

FX100
FLOPS /
Time [s]
PEAK
[%]
187

1.64

181

1.69

128

2.38

126

2.44

100

2.95

コンパイラによる SIMD 化促進のため指示行を挿入(ウ) 複素数演算のチューニング 2a)+高コストな多重ループで最内側ループ長

前項とは別のケースとして，FX10 上にてノーが長くなるようにループ入れ替え

ドあたりの計算負荷を減らしたケースで複素数b)+realloc 呼び出し回数の削減

演算のチューニングを行った．使用する
5

学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月

AdvMag のバージョンは 1.6.1 であり，モジュー

2000.

ルはハイブリッド並列版を用いる．複素数の実[2] 宮村倫司，野口裕久，塩谷隆二，吉村忍，

装は C99&Macros とし，コンパイラオプション

矢川元基，階層型領域分割法による超並として-Kfast を用いてコンパイルしたものを基

列弾塑性有限要素解析，日本機械学会論本とする．モデルとして 100 万自由度のケーキ

文 集 (A 編 ), 65-634, pp.1201-1208 ，モデルを用い，部分領域数を 8,800 とする．計1999.

測には 8 ノードを用いる．

6.

本項での検討として，以下を行う．
a)

今年度の進捗状況と今後の展望

最内側ループを行で回していた 2 次元配列を

ほぼ計画通りに進めることができた．ピーク

扱うループで，最内側ループを列で回すよう性能比もチューニング前は FX10 で約 0.8%であにループ入れ替えったが，最終的に 7.41%まで向上させることがでb)a) +realloc 呼び出し回数の削減

きた．ただし FX100 では半分程度に落ちてしまc)b) +コンパイラによる SIMD 化促進のためのい，課題が残った．これは SIMD 化が不十分で，指示行を挿入FX10 から 2 倍になった SIMD 演算の幅に対応で

c) +コンパイラオプションによる最適化を試きていていないためと思われる．今後はさらに

行：-Kfast,ocl,restp=arg,preex

SIMD 化を進める検討を行うとともに，より高並

表 3 にチューニング結果を示す．ノードあた

列な環境での性能向上を検討していく．

d)

りの計算負荷が小さいにも関わらず前項よりも計算時間が長くなっているのは，テストした反7.

研究成果リスト

復回数が異なるためである．ピーク性能比は(1) 学術論文

1.24 倍程度にとどまったが，ノードあたりの計[1] Amane

TAKEI,

Kohei

MUROTANI,

算負荷が小さいにも関わらず前項よりもピーク

Shin-ichiro SUGIMOTO, Masao OGINO

性能比が高く，最終的に 7.41%まで向上するこandとができた．チューニング開始前は約 0.8%であElectromagnetic Field Simulation using

ったことから，一連のチューニングによって 9

Numerical Human Body Models, IEEE

倍以上の効率化を達成できた．

Transactions on Magnetics, in press.

Hiroshi

KAWAI,

High-Accuracy(2) 国際会議プロシーディングス
表 3. チューニング結果 2

asis
a)
b)
c)
d)

Time [s](ratio to asis)
850.63(1.00)
725.48(1.17)
703.35(1.21)
690.18(1.23)
683.17(1.25)(3) 国際会議発表[1] Shin-ichiro

FLOPS / PEAK [%]Daisuke

TAGAMI, Masao OGINO, Amane TAKEI

5.95

and Hiroshi KANAYAMA, Improvement

6.96

of Convergence in Time-Harmonic Eddy

7.20

Current Analysis by Hierarchical Domain

7.31

Decomposition

Method,

7.41

International

Conference

The

20th

on

the

Computation of Electromagnetic Fields(COMPUMAG2015), USB, 2015.

参考文献[2] Amane[1] 金山寛，岩波講座 現代工学の基礎 ＜空間系 IV＞SUGIMOTO,TAKEI,

Kohei

MUROTANI,

Shin-ichiro SUGIMOTO, Masao OGINO

計算電磁気学, 岩波書店，

and6

Hiroshi

KAWAI,

High-accuracy

学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月

electromagnetic field simulation using

法ソルバの高度化・高速化の取り組み，日numerical human body models, The 20th

本 機 械 学 会 第 28 回 計 算 力 学 講 演 会

International(CMD2015)，CD-ROM，2015.

Conference

on

the

[7] 金山寛，荻野正雄，杉本振一郎，鄭宏杰，Computation of Electromagnetic Fields(COMPUMAG2015), USB, 2015.

淀薫，静磁場領域分割インターフェイス問[3] Hiroshi KANAYAMA, Masao OGINO and

題の前処理付き共役勾配法

Shin-ichiro SUGIMOTO, Preconditioners[8] 武居周，領域分割法に基づく高性能並列電in Domain Decomposition Method for

磁界解析手法，日本機械学会 第 28 回計算

Magnetostatic

力学講演会(CMD2015)，CD-ROM，2015.

International

Problems,

The

Conference

6thon(5) その他（特許，プレス発表，著書等）

Computational Methods, 2015.
[4] Shin-ichiro SUGIMOTO, Implementation
and Acceleration of the HDDM for theElectromagnetic

Field

Problemreport),

34th

JSST

The

(2nd

AnnualConference, International Conference on
Simulation Technology, online, 2015.
(4) 国内会議発表[1] 金山寛，荻野正雄，杉本振一郎，鄭宏杰，
大規模静磁場インターフェイス問題用前処理の開発，第 27 回「電磁力関連のダイナミクス」シンポジウム講演論文集，pp.451-454，
2015.
[2] 武居周，杉本振一郎，ボクセルメッシュモ
デルを用いる並列高周波電磁界解析，第 27回「電磁力関連のダイナミクス」シンポジ
ウム講演論文集，pp.455-458，2015.
[3] 杉本振一郎，領域分割法における時間調和
渦電流解析の収束性改善，第 27 回「電磁力関連のダイナミクス」シンポジウム講演論
文集，pp.459-460，2015.
[4] 武居周，杉本振一郎，ボクセルメッシュモ
デルを用いる並列 full-wave 電磁界解析，計算工学講演会論文集，Vol.20，CD-ROM，
2015.
[5] 杉本振一郎，田上大助，荻野正雄，武居周，
金山寛，階層型領域分割法における時間調
和渦電流解析の収束性改善(第 2 報)，計算工学講演会論文集，Vol.20，CD-ROM，2015.
[6] 杉本振一郎，電磁場解析向け並列有限要素7