EEM-FDM 電磁界シミュレータ〇差分法を用いた電磁界シミュレータ時間領域差分法 (FDTD 法 ) : 高周波用 ( 計算対象 > 波長 ) 周波数領域差分法 : 低周波用 ( 計算対象 < 波長 ) 〇用途準静電界 (10 3 Hz) から光 (10 15 Hz) まであらゆる電磁界現象

Size: px

Start display at page:

Download "EEM-FDM 電磁界シミュレータ〇差分法を用いた電磁界シミュレータ時間領域差分法 (FDTD 法 ) : 高周波用 ( 計算対象 > 波長 ) 周波数領域差分法 : 低周波用 ( 計算対象 < 波長 ) 〇用途準静電界 (10 3 Hz) から光 (10 15 Hz) まであらゆる電磁界現象"

しおりかんけ
5 years ago
Views:

1 EEM-FDM のご紹介 ~GPU 高速計算を中心に ~ 株式会社 EEM 2017/03/14

2 EEM-FDM 電磁界シミュレータ〇差分法を用いた電磁界シミュレータ時間領域差分法 (FDTD 法 ) : 高周波用 ( 計算対象 > 波長 ) 周波数領域差分法 : 低周波用 ( 計算対象 < 波長 ) 〇用途準静電界 (10 3 Hz) から光 (10 15 Hz) まであらゆる電磁界現象に対応アンテナ EMC 導波管マイクロストリップ線路マイクロ波回路電波伝播電波吸収体電波暗室メタマテリアル〇高速化技術 FDTD 法は並列化に適しているので CPU/GPU( クラスタ含む ) により大規模計算が可能〇動作環境 WindowsPC からスパコンまで NVIDIA のグラフィックスボードを用いた超高速計算

3 FDTD 法 (Finite Difference Time Domain 時間領域差分法 ) Maxwell の方程式 : E = B t H = D t +J D=ϵ E B=μ H J =σ E E : 電界 D : 電束密度 H : 磁界 B : 磁束密度 J : 電流 ϵ: 誘電率 μ : 透磁率 σ : 導電率 D,B,J を消去 E = μ H t H=ϵ E t +σ E E,H: 未知数 ε,μ,σ: 既知数宇野亨 FDTD 法による電磁界およびアンテナ解析コロナ社,

4 時間領域で離散化 (LeapFrog : 蛙跳びアルゴリズム ) Z H n+ 1 2 =d 1 Z H n 12 d 2 (c Δ t) E n E n+1 =c 1 E n +c 2 (c Δ t) Z H n+ 1 2 c 1, c 2, d 1, d 2 : 電気定数で決まる定数 ( 場所の関数無次元 ) Z = 120π [Ω] c : 光速 H n 1 2 Δt n+ 1 2 E n Δt n+1 時間 4

5 空間領域で離散化 (X 成分左辺が最新項 ) 並列化に極めて適している!! d 2 ( i, j+ 1 2, k ){E z +c 2( i+ 1 2, j, k ) {ZH n+ 1 ( 2 z ZH x n+ 1 2( i, j+ 1 2, k ) = d 1( i, j+ 1 2, k ) ZH x n( i, j+1, k + 1 ) 2 E n ( i, j, k+ 1 ) z 2 E x n+1( i + 1 2, j, k ) = c 1( i + 1 2, j, k ) n( E i+ 1 x 2, j, k ) i + 1 2, j+ 1 2, k ) ZH z Z Δ y j+ 1 2 c Δ t Δ y j c Δ t n+ 1 2 ( i+ 1 2, j 1 2, k ) n 1 2( i, j+ 1 2,k ) n( E i, j + 1 y 2,k +1 ) E n ( i, j+ 1 y 2, k ) } Δ z 1 k + 2 c Δ t 2( 1 n+ ZH i+ 1 y 2, j,k + 1 ) 2( 2 ZH n+ 1 i+ 1 y 2, j,k 1 ) 2 } Δ z k c Δ t Hy Ez Hx Ey Y Yee 格子電界 E: 稜線の中心の接線方向磁界 H: 面の中心の法線方向 Ex Hz 5 X

6 吸収境界条件計算領域吸収境界条件 : 計算領域境界での電磁波の反射をなくす Mur 一次 : 1 層計算対象 PML : 5~10 層 ( メモリーと計算時間が少し増えるが精度が向上する ) 6

7 フーリエ変換により時間領域から周波数領域に変換する E (r,ω) = E(r, t)e j ωt dt 0 H (r,ω) = 0 入力インピーダンス H (r,t)e j ω t dt Z in (ω)= V in(ω) I in (ω) (V は E から I は H から計算する ) 遠方界 j k r j k e E { θ (r,θ,ϕ) = ϕ } 4 π r F { θ (θ,ϕ) ϕ } F { θ ϕ} (θ,ϕ) = ZN { θ ϕ} (θ,ϕ)±l { ϕ (θ,ϕ) θ} N (θ,ϕ) = S L(θ,ϕ) = S ^n H (r)exp( j k r ^r )ds ^n E(r)exp( j k r ^r)ds (S: 計算領域の境界面 ) 7

8 EEM-FDM の新機能新機能 PC PC PC PC メモリーメモリーメモリーメモリー CPU CPU CPU CPU CPU CPU CPU コア OpenMP OpenMP LAN MPI PC PC PC PC メモリーメモリーメモリーメモリーメモリーメモリーメモリー GPU GPU GPU GPU GPU GPU GPU コア CUDA CUDA+MPI LAN CUDA+MPI 8

9 ノードプロセススレッドノードノードスレッドプロセスプロセスプロセスプロセススレッドネットワークノード : CPU とメモリー ( と GPU) を持っている : PC 一台と同じプロセス : 独立したメモリー空間を持っている独立したプログラムスレッド : プロセスが管理する一連の計算処理 ( コアとほぼ同義 ) クラスタは複数のノードから成りネットワークでつながっているノードは複数のプロセスを起動する ( タスクマネージャーに表示される ) プロセスは複数のスレッドを実行する逆に複数のプロセスにまたがるスレッドはない複数のノードにまたがるプロセスはない複数のノードとプロセスが起動できるのは MPI のみ 9

10 並列化プログラミング新機能 OpenMP CPU 向け共有メモリーモデル領域分割が不要プログラミングが容易 MPI CPU クラスタ向け分散メモリーモデル領域分割が必要プログラミングが複雑 CUDA GPU 向け共有メモリーモデル領域分割が不要プログラミングが複雑 CUDA+MPI マルチ GPU と GPU クラスタ向け分散メモリーモデル領域分割が必要 CUDA と MPI が適切に実装されていればプログラミングは比較的容易 MPI : Message Passing Interface 分散メモリー向けのプログラミング仕様共有メモリー環境でも動作しパフォーマンスも落ちない MPI で実装すればすべての環境で最高のパフォーマンスが得られる ( ただしクラスタでは高速なネットワーク環境が必要 ) 開発環境はすべて無償 10

11 領域分割と MPI 〇 FDTD 法はアルゴリズム上領域分割が必要になる〇 MPI の通信 : タイムステップ毎に不足成分を隣のプロセスから通信によって取得する不足成分 : 領域境界の半セル外側の Hy,Hz 重複成分 : 領域境界上の Ey,Ez,Hx 領域境界通信方向 Ey Ez Hy Hz Hx Y プロセス P プロセス P+1 X 11

12 EEM-FDM の必要リソースメモリーメモリー = (30 / プロセス数 + 24 * 周波数数 ) * Nx * Ny * Nz [ バイト ] 例えばプロセス数 =1, 周波数数 =1 のとき Nx=Ny=Nz=100 のとき 54MB Nx=Ny=Nz=1000 のとき 54GB ファイルサイズ入力ファイルサイズ : 一般に小さい出力ファイルサイズ = ( * 周波数数 ) * Nx * Ny * Nz [ バイト ] 使用メモリーとほぼ同じ SSD 推奨十分な空き容量が必要計算時間計算時間 Nx * Ny * Nz * タイムステップ数 / プロセス数 12

13 グラフィックスボード ( ビデオカード GPU) NVIDIA のグラフィックスボードが CUDA に対応現在の世代 :GTX 10XX (Pascal 世代 = Compute Capability 6.X) 13

14 FOCUS スパコン ( 公益財団法人計算科学振興財団 ) 14

15 FOCUS スパコンの構成システム CPU コア数 ( ノード当たり ) メモリー ( ノード当たり ) ノード数ノード間ネットワーク使用料 (1 ノード 1 時間 ) B Xeon E GB 2 D Xeon E5-2670v GB 80 E Xeon E5-2670v GB 48 F Xeon E5-2698v GB 12 H Xeon D GB 136 InfiniBand 40Gbps InfiniBand 56Gbps InfiniBand 56Gbps InfiniBand 56Gbps Ethernet 10Gbps x 円 300 円 300 円 500 円 100 円アカウント料年間 1 万円が必要使用料はジョブ従量制類似のサービス :Amazon Web Services (EC2) : 15

16 計算時間の一例セル数 =300 X 300 X 300 タイムステップ数 = 計算時間 [ 秒 ] 10 Xeon D-1541 Xeon E5-2670v2 i7-4770k GTX プロセス数,GPU 数 16

17 クラスタは何台 ( 何プロセス ) まで速くなるか通信時間が計算時間と同じオーダーになると速度比が低下し始める 1 タイムステップ当たりの計算時間 = (Nx * Ny * Nz) / (300 * 300 * 300) * (500 / 500) / Np (i7-4770k の場合 ) = 37 * 10-9 * (Nx * Ny * Nz / Np) [ 秒 ] (Np: プロセス数 ) 通信時間 = 遅延時間 + 通信量 / 通信速度 = (2 * Ny * Nz * 4 * M) / (bps / 8) [ 秒 ] ( 遅延時間 ~10μsは無視できる ) 通信時間 << 計算時間となるためには Np << Nx / 7 * 通信速度 [Gbps] (M: 両側双方向通信のため M=4 としている ) 結論 : クラスタの有効台数は Nx と通信速度に比例する例えば 1Gbps, Nx=300 のとき Np << 40 プロセスまで (~Beowulf クラス ) Nx に比例する理由は X 方向に領域分割しているため計算領域が細長いときは X 方向を長手にすると有利上の評価式が成り立つ条件 : Nx >> Np 領域境界での重複計算吸収境界条件 17

18 クラスタの使い方準備作業 1. 計算に使用するすべてのコンピュータをネットワークで接続します ( 高速なネットワーク推奨 ) 2. すべてのコンピュータに同じアカウントとパスワードを設定します ( パスワードなしは不可 ) 3. すべてのコンピュータの同じ絶対パスに 4 つのファイル (fdm2_mpi.exe, fdm2g_mpi.exe, msmpi.dll, smpd.exe) をコピーします計算 1. [ ホスト名 ] と [ プロセス数 ] を適当に設定します CPU クラスタでは [ プロセス数 ] に物理コア数と論理コア数の間の数値を入力してください GPU クラスタでは [ プロセス数 ] にグラフィックスボードの数を入力してください 2. すべてのコンピュータでコマンドプロンプトで smpd -p 8677 を実行します (8677 は MPI のポート番号 ( 固定 )) 3. [ 計算 ] ボタンをクリックすると計算を開始します 18

19 データ作成ライブラリー以下のようなケースではデータ作成ライブラリーが便利 (1) データの数が多くかつ規則的なとき (2) 一つの座標を変えると同時にたくさんの変更が必要になるとき (3) パラメーターを変えながら繰り返し計算したいとき右に使用例を示す ( ダイポールアンテナ ) 用意された関数を呼び出してデータを作成する C 言語の初歩的な知識が必要 #include "fdm_datalib.h" int main() { // (1) initialize fdm_init(); // (2) title fdm_title("sample1"); // (3) domain fdm_domain(0); // (4) mesh fdm_xsection(2, -50e-3, +50e-3); fdm_xdivision(1, 20); fdm_ysection(2, -50e-3, +50e-3); fdm_ydivision(1, 20); fdm_zsection(4, -75e-3, -25e-3, +25e-3, +75e-3); fdm_zdivision(3, 10, 11, 10); // (5) material fdm_material(2.0, 0.0, 1.0, 0.0); // (6) geometry fdm_unit6(1, 1, 0e-3, 0e-3, -25e-3, 0e-3, 0e-3, 25e-3); // dipole // (7) incidence fdm_feed(3, 0e-3, 0e-3, 0e-3, 1, 0); // (8) observation point // (9) frequency fdm_freq(3e9, 3e9, 0); fdm_freq2(2e9, 4e9, 20); // (10) solver fdm_iteration(1e-3, 2000, 100); // (11) output fdm_outdat("sample1.fdm"); 19 } return 0;

20 計算例 20

21 パッチアンテナ (1/2) パッチアンテナ計算モデル ( ミリ波用 ) パッチ =2mm 基板厚さ =0.4mm 基板比誘電率 =4 反射損失周波数特性 (30~35GHz, Z0=50Ω) 21

22 パッチアンテナ (2/2) アンテナ単体金属板の障害物 300mm 400mm 電界分布 (33GHz) -20dB +10dB 放射パターン (33GHz) セル数 = 660 x 660 x 455 = 2.0 億タイムステップ数 =7500 計算時間 =14 分 (@2GPU) 使用メモリー =6.2GB x 2 GPU : NVIDIA GeForce GTX 1070 (8GB GDDR5) X 2 ( 以下同じ ) 22

23 車載アンテナアンテナ ( 垂直モノポール右端 ) 後前電界分布 (6GHz 中心面 ) セル数 = 440 x 1100 x 404 = 2.0 億タイムステップ数 = 2300 計算時間 = 5 分 (@2GPU) 使用メモリー = 6.2GB x 2 左右放射パターン (6GHz) 23

24 新幹線セル数 = 2141 x 322 x 217 = 1.5 億タイムステップ数 = 6700 計算時間 = 10 分 (@2GPU) 使用メモリー = 4.0GB x 2 アクセスポイント電界分布 (2.45GHz, 床上 80cm) 24

25 航空機セル数 = 400 x 2560 x 135 = 1.4 億タイムステップ数 = 5000 計算時間 = 9 分 (@2GPU) 使用メモリー = 3.8GB x 2 AP AP 電界分布 (1.5GHz) 25

5GHz,Z=6mm) BBE(Backfire-Broadside-Endfire) アンテナ : 周波数でビームの向きが変わるアンテナ C.Caloz and T.

26 メタマテリアル C L stub inductor interdigital capacitor メタマテリアルの等価回路吸い込み口湧き出し口マイクロストリップ線路 via 進行方向 Ez の位相分布 (2.5GHz) 左手系 S11,S21 の周波数特性 (0~6GHz) Ez の位相分布 (2.5GHz,Z=6mm) BBE(Backfire-Broadside-Endfire) アンテナ : 周波数でビームの向きが変わるアンテナ C.Caloz and T.Itoh, Electromagnetic Metamaterials: Transmission line theory and Microwave Applications, Wiley Interscience,

27 ダブルリッジホーンアンテナ E 面電界分布 (3GHz) E 面放射パターン (3GHz) S11 周波数特性 (2~8GHz) H 面電界分布 (3GHz) H 面放射パターン (3GHz) 27

28 電波暗室 (1/2) 電波暗室の大きさ = 10m x 5m x 5m セル数 = 1000 x 500 x 500 = 2.5 億 ( セルサイズ =1cm) 使用メモリー = 6.6GB x 2 (1) 周波数 = 300MHz タイムステップ数 = 2900 計算時間 = 8 分 (@2GPU) (2) 周波数 = 1GHz タイムステップ数 = 1900 計算時間 = 6 分 (@2GPU) 28

29 電波暗室 (2/2) 電波暗室内の電界分布 300MHz 垂直面 300MHz 水平面床上 1.5m 1GHz 垂直面 1GHz 水平面床上 1.5m 29

s ss s ss = ε = = s ss s (3) と表される s の要素における s s = κ = κ, =,, (4) jωε jω s は複素比誘電率に相当する物理量でありここで PML 媒質定数を次のように定義するすなわち κξ をPML 媒質の等価比誘電率 ξ をPML 媒質の

s ss s ss = ε = = s ss s (3) と表される s の要素における s s = κ = κ, =,, (4) jωε jω s は複素比誘電率に相当する物理量でありここで PML 媒質定数を次のように定義するすなわち κξ をPML 媒質の等価比誘電率 ξ をPML 媒質の FDTD 解析法 (Matlab 版 2 次元 PML) プログラム解説 v2.11 1. 概要 FDTD 解析における吸収境界である完全整合層 (Perfectl Matched Laer, PML) の定式化とプログラミングを2 次元 TE 波について解説する PMLは異方性の損失をもつ仮想的な物質であり侵入して来る電磁波を逃さず吸収する通常の物質と接する界面でインピーダンスが整合しており