2014/4/25
先駆的科学計算に関するフォーラム2014(Forum on Advanced Scientific Computing 2014)〜先端的計算科学研究プロジェクト成果報告〜

並列有限要素法に基づく人体内部の高周波電磁界の大規模シミュレーション

武居周（苫小牧高専 電気電子工学科）

目次
 イントロ，目的
 基礎式，有限要素定式化，方程式の性質，領域分割法
 解析：数値人体モデルを用いた人体内部の電波伝搬解

析 性能評価 まとめ

イントロ：電磁界シミュレーション


国内における研究例

アンテナ設計，電気機器設計，電磁環境問題の評価などに適用

パラボラアンテナのシミュレーション（株式会社フォトン）

モーターのシミュレーション（北海道大学）鉄道車内において携帯電話を使用した電磁界解析（苫小牧高専：本研究による成果のひとつ）

イントロ：電磁界シミュレーション

ETCシステムやキーレスエントリーシステムなど
自動車車内および周辺において本研究において着目する対象
通勤電車などの公共の場における携帯電話の使用

電磁界温熱治療予測技術向上に向けた人体内部の電波伝搬解析（出典：山本ビニターHP: http://www.vinita.co.jp/）

生活環境中全体にわたる微弱な高周波電磁界の評価手法，電磁界を利用する医療技術等への適用の要求の高まりに対応したい

本研究の目的
 人体内部などの大規模かつ複雑なモデルによる高精度な計算が可能な高周波電磁界解析手法の検討
本研究のねらい:
階層型領域分割法を用いる高周波電磁界の大規模
解
析手法の研究・開発推進，および環境モデルや人体モデルをはじめとする実用的な数値モデルを用いる

基礎式
 変位電流を含むMaxwell方程式より導かれる電界Eを未知数とするベクトル波動方程式を基礎式とする
curl (m -1 curlE) - w 2e E - jws E = jw J in W

E  n  0 on  ( curlE) ´ n - jw e0 / m0 ( E ´ n) ´ n = 0on ¶WABCJ σEn∂
解くべき領域

Notice:E [V/m]: 電界（未知数）, E [V/m] : 既知節点値

ω=2πf [rad/s]: 角周波数（ f [Hz]：周波数）

ε0 [F/m]: 真空の誘電率, εr : 比誘電率μ0 [H/m]: 真空の透磁率, μr : 比透磁率

J [A/ m2 ]: 電流密度，σ[S/m]：導電率

E法による有限要素定式化
 電界Eを未知数とするベクトル波動方程式の弱形式を導き，辺要素によって以下の有限要素方程式を得る-1*2**mcurlE×
curlEdvw
eE×
EdvjwsE×Eòòò W
òòò W h hòòò W h hdv
hh+ jw e0 m0 òò¶WABC(Eh ´ n)×(E*h ´ n)dS = jw òòò J h ×E*hdv,W

Notice:n∂自由度

Eh: Eの有限要素近似

Eh *: Eh*×n=0となる任意の試験関数

Jh: Jの有限要素近似

領域分割法による並列化

領域分割法(Domain Decomposition Method: DDM)
 解析領域（要素分割された有限要素モデル）を空間的にサブ領域（部分領域）に分割することで，演算を部分領域の独立領域と部分領域間の袖領域に分けて実行する手法

領域分割は要素ベースで行う部分領域で要素の重なりはない
有限要素モデル
有限要素辺自由度

領域分割モデル
領域独立辺

袖領域辺

DDMアルゴリズム
 解くべき行列式 ku = f の分割


内部(I), 内部境界(B)で自由度を並び替えた式
é(1)ê K IIêê0êêê RB(1) K II( N )Têëù
é
ù
fI(1)
úé u(1) ù êúIêúúêúêúúêú
(N )=K IB( N ) RB( N )fê
úIú u( N )êú
Núê I ú ê Nú(i)T(i) ( N )(i) (i)êúå RB K BB RB úë uB û ê å RB fB úúûêë i=1úû
i=1
K IB(1) RB(1)0K II( N )RB( N ) K IB( N )T

部分領域内部と領域間境界に行列方程式を分割(1) 領域間境界上自由度に関する式（インターフェース問題）
1

1

M SuB  M g ,N

S = åR(i)T
B

i=1
N

T(i )
(i )
S (i )  K BB
 K IB(i ) ( K II(i ) ) 1 K IB

S R ,(i)(i)Bg   RB(i ) f B( i )  K IB(i ) ( K II( i ) ) 1 f I( i )i 1(2) 部分領域内部自由度に関する式

K II(i ) u I(i )  f I(i )  K IB(i ) u B(i ) ,i  1,..., NT(1) を解き，求めたuBを用いて(2)を解き，全体の解uを求める

DDMアルゴリズム
893/ : ' $; <$8<

計算の流れ

J %5$?%/ - ' 1: ') *+&%2 3,/ $$) *+&%2 3,/ $7) *+&%2 3,/ $6) *+&%2 3,/ $$) *+&%2 3,/ $7) *+&%2 3,/ $6) *+&%2 3,/ $$) *+&%2 3,/ $7) *+&%2 3,/ $64315$64315$64315$6
M' ?' ,- ' $$M' . *(5.

M' 3&$> 353

8%/ - ' 1: '"J >A/ ' $2 %1' $3/ 3(O. ,. $%0$) *+&%2 3,/ . $4315$74315$74315$74315$$

4315$$

4315$$> %2 3,/ $> ' ?%2 C%. ' 1=$
H> I "J KL M"N# "K@)) *+&%2 3,/ $) %(- ' 1=$!" ! #$> ,1' ?5$. %(- ' 1@/ 5' 103?' $) %(- ' 1=$
8A8B$# ' 59%&

! ,' 131?9,?3($" %2 3,/ $" ' ?%2 C%. ,D%/ $# ' 59%&$EF > > # G

NICT数値人体モデルを用いた人体内部の電波伝搬解析


癌温熱治療装置の治療計画立案に役立てることを目的とす
る，NICT数値人体モデルを用いた人体内部の電波伝搬解析



これらのモデルをNICT（情報通信研究機構）が主導し開発
人体（組織・臓器）の形状を微小のボクセルで表現
日本人成人男女の平均身長と平均体重に合致したボランティアのMRIデータに基づき、2mmの空間分解能、51の組織・臓器を有する
無償（非営利目的）で公開
ボクセル数，四面体分割時の要素数•辺数
右：男性モデル(TARO)，左：女性モデル（HANAKO）
 320x160x866(TARO) = 4,400万 voxels
出典：情報通信研究機構(NICT) Homepage
http://emc.nict.go.jp/bio/model/model01_1.html
 320x160x804(HANAKO) = 4,100万 voxels
 全てを四面体にすると、2.2億 tetras
 辺の総数は、2.6億（267,437,186）ボクセルを四面体に分割

ボクセルを四面体に分割

四面体に分割ボクセルデータからHDDM
有限要素メッシュを生成するツール
メッシュ生成アルゴリズム
1. 数値人体データ（ボクセル）の読み込み，双対グラフの作成
2. part分割後に，各partのボクセルを四面体分割
3. 各partにおいてsubdomain分割
Voxel
双対グラフ
METIS
四面体分割

四面体分割

四面体分割

四面体分割

METIS

METIS

METIS

METIS

SolverSolverSolverSolver

数値例
 人体モデルに適用する生体物性パラメータ• （独）理化学研究所 生体力学シミュレーション身体組織物性値データベース
http://cfd-duo.riken.go.jp/cbms-mp/index.htm※1で示したものは女性(Hanako)用の生体物性パラメータであり、数値モデルに存在しない
ID4, 37, 42とあわせて対象から除外してある※2で示した血液および臓器の内容物に関して、本データベースに無いため、生理食塩水と同等程度の導電率および比誘電率とした

数値例
 人体モデルに適用する生体物性パラメータ• （独）理化学研究所 生体力学シミュレーション身体組織物性値データベース※1で示したものは女性(Hanako)用の生体物性パラメータであり、数値モデルに存在しない
ID4, 37, 42とあわせて対象から除外してある※2で示した血液および臓器の内容物に関して、本データベースに無いため、生理食塩水と同等程度の導電率および比誘電率とした

HPC環境における性能評価
 数値モデル NICT数値人体モデル（TARO, 全身）→Mesh 2: 4mmピッチで分割したモデル
2,771万要素
3,190万複素自由度→Mesh 4: 2mmピッチで分割したモデル
2.2億要素2.6億複素自由度
右：男性モデル(TARO)
左：女性モデル（HANAKO）
出典：情報通信研究機構(NICT) Homepage
http://emc.nict.go.jp/bio/model/model01_1.html

HPC環境における性能評価
 HPC環境

九州大学情報基盤センターFX10
 Fujitsu PRIMEHPC FX10 x 8 racks
 768ノード，181.6 PFlops
 6次元メッシュ/トーラスTofuインターコネクト
 プロセッサ:SPARC64 Ix-fx
 ノード：


16 cores of 1.848 Ghz
and 32 GB memory

HPC環境における性能評価


FX10における，実用モデル（NICT数値人体モデル）による
DDM-COCGの性能評価Ⅰ
Calculation Time and Average Memory Requirement
Model

Mesh 2Mesh 4ComplexDOFs

27.21 M

260 M

Memorysize / node(MByte)

160

160

Num. ofCores

Frequency(MHz)

Iteration
countsCalculation
time (hour)

12,8752.7

320,81,641

1.4

20node70

2,199

1.9300

9,684

8.5

19,0287.7

320,84,817

4.520node

70

7,411

6.3

300

18,30615.6

解析結果の可視化例

人体内部の電界E[V/m]

人体内部の磁界H[A/m]
Mesh4の解析結果の可視化(300[MHz])

HPC環境における性能評価
FX10における，実用モデル（NICT数値人体モデル）による
DDM-COCGの性能評価Ⅱ
Strong Scaling of the HDD Mesh of 4mm (Mesh 2)
Model

Num. ofparts(=Cores)

HDDM
loop (sec)

Speedup

Ideal speedup

Mesh 2 -1

320

227

1.000

1-2

640

144

1.576

2-3

1280

102

2.225

4-4

2560

66

3.439

8-55120425.40516

ノード数2倍で1.575倍の高速化

HPC環境における性能評価
FX10における，実用モデル（NICT数値人体モデル）による
DDM-COCGの性能評価Ⅱ
Strong Scaling of the HDD Mesh of 2mm (Mesh 4)
ModelNum. ofparts(=Cores)

HDDM
loop (sec)

Speedup

Ideal speedup

Mesh4 -1

320

2195

1.000

1-2

640

1293

1.698

2-3

1280

866

2.535

4-4

2560

548

4.005

8-5

5120

365

6.014

16

ノード数2倍で1.698倍の高速化

Mesh 2よりも問題規模が大きいMesh 4の方が、スピードアップ率が大きい。

PCクラスタ（Core i7, 80core）
を使用した計算例
Model

Num. oftet. elements

Memorysize[G Byte]PE

Num. of
Subdomain
per a Part

Mesh 227.7 M

0.180

4,000

Mesh 4

220 M

0.180

23,000Analysis
frequency[MHz]

Iteration
countsCalc.
time[h]1870300
1870300

2,875
1,641
2,199
9,684
9,028
4,817
7,411
18,306

1.5
0.8
1.1
4.8
20.4
11.0
16.8
41.4

• もっとも計算時間を要するMesh 4の300MHzにおいて41.4[h]で計算が完了• 数億要素規模では，PCクラスタの方が有利

HPC環境向けチューニングが必要（実施中）複素数実装の検討
 複素数を扱うHDDMでの新たな実装方法を検討
 データの表現方法• 構造体(struct Complex { double re; double im; })• double complex型(C99) 四則演算の実装方法

• 関数
• マクロ
 従来の実装と比べて，4倍の高速化

データの表現方法
 高周波電磁界解析コードにおける従来の表現方法
 構造体“Complex”
struct Complex {double re ; /* 実部 */double im ; /* 虚部 */};

 新たな表現方法
 C99規格(ISO/IEC 9899:1999)で定められている“double complex”型

四則演算の実装方法“Structure & functions”
void ComplexPlus( Complex *x, Complex a, Complex b ){(*x).re = a.re + b.re ;(*x).im = a.im + b.im ;} /* addition */void ComplexMulti( Complex *x, Complex a, Complex b ){(*x).re = a.re*b.re - a.im*b.im ;(*x).im = a.im*b.re + a.re*b.im ;} /* multiplication */“Structure & macros”#define ComplexPlus( __x, __a, __b ) { ¥(*__x).re = __a.re + __b.re ; ¥(*__x).im = __a.im + __b.im ; }#define ComplexMulti( __x, __a, __b ) { ¥double __re = __a.re*__b.re - __a.im*__b.im ; ¥(*__x).im = __a.im*__b.re + __a.re*__b.im ; ¥(*__x).re = __re ; }“double complex & macros”#define ComplexPlus( __x, __a, __b ) ¥{ *__x = __a + __b ; }#define ComplexMulti( __x, __a, __b ) ¥{ *__x = __a * __b ; }

計算条件
 計測
 計算時間[s]，ピーク性能比[%] 計測はプログラム全体で行う• データの入力，解析，結果出力など
 比較
 old: “structure & functions”
 new: “structure & macros”

 C99: “double complex & macros”

HPC環境向けチューニング
old: “structure & functions”（構造体＋四則演算等の関数）
new: “structure & macros”（構造体＋関数のインライン展開）
C99: “double complex & macros”（double complex型の利用）
ピーク性能比 [%]
Nodes

old

new

C99

12

1.191

5.010

24

1.185

4896

Speed-up
new

C99

7.244

4.21

6.08

4.957

7.169

4.18

6.05

1.171

4.818

6.868

4.12

5.87

1.148

4.599

6.619

4.01

5.76

インライン展開によって，ピーク性能比が改善 double complex型が構造体に比べて45%高いまとめ 高周波電磁界コード 大規模解析手法の開発を継続的に実施することが必要

 基礎式，有限要素定式化，領域分割法
 E法によるベクトル波動方程式

 解析例：人体内の電波伝搬解析
 2.2億要素の全身モデルの解析に成功

 HPC環境における性能改善が課題まとめ 今後の予定（展望を含む）
 HPC環境向けにコードを改良（実施中）• 本解析コードを利用し，少なくとも100億要素程度の高周波電磁界解析が可能なコード開発を目指す• 本解析コードのチューニングにより，ピーク性能比をさらにあげてゆく 人体シミュレーション機能の，癌温熱治療シミュレー

タ構築に向けた高機能化開発（実施中）• 熱効果を表現するための熱問題との連成機能（開発・検証）• ボクセル形状に起因するノイズ除去のためのメッシュスムージング機能（開発中）• 呼吸運動を表現するための構造問題との連成機能（検討中）• 血流冷却モデルの導入（検討中）Thank youご清聴ありがとうございました。

Furano town, Hokkaido , Japan