JP2008052530A

JP2008052530A - プラズマ粒子シミュレーション計算方法

Info

Publication number: JP2008052530A
Application number: JP2006228728A
Authority: JP
Inventors: Katsumi Ichifuji; 克己一藤; Iwao Ryu; 岩劉
Original assignee: HAMAMATSU METRIX KK
Current assignee: HAMAMATSU METRIX KK
Priority date: 2006-08-25
Filing date: 2006-08-25
Publication date: 2008-03-06

Abstract

【課題】プラズマ粒子シミュレーションVillasenor-Buneman法の８０％を占める計算ボトルネック部を高速に計算することができるプラズマ物理の計算機シミュレーション計算方法を提供する。
【解決手段】空間メッシュの格子立方体と同じ体積の粒子を想定し、この粒子に一定密度の電荷を帯びさせ電磁場の空間に配置することでプラズマの挙動を模擬する、Villasenor-Buneman（ヴィラセノウ−ブーネマン）法を汎用計算機を用いたプラズマ物理の計算機シミュレーションにおいて、Newton-lorentz（ニュートン−ローレンツ）の運動方程式（式（１））の計算を、電流密度の計算（式（２））、Maxwell（マックスウェル）の方程式の計算（式（３））、粒子に働く力の計算（式（４））、と段階的に計算をするに際し、前記式（２）の計算を実行させるための計算ユニットを前記汎用計算機内に実装して、プラズマ粒子の挙動を模擬する。
【選択図】図８

Description

本発明は、プラズマ物理の計算機シミュレーションに関し、さらに詳しくは、プラズマ粒子シミュレーション計算方法に関する。

プラズマ物理の計算機シミュレーションを実施する場合、一般に使われる粒子シミュレーション法は、ＰＩＣ（Particle In Cell）と呼ばれ、粒子間の相互作用を直接解くのではなく、Newton-Lorentz運動方程式を用い、粒子の座標、速度などから空間格子状の電荷・電流を求めPoisson方程式等で空間電場・磁場分布を計算しMaxwell方程式で粒子に加わる力を計算するという手法を用いている。
しかし、微分形式からみた場合、このような操作は差分形式では極めて複雑な問題となるので、Villasenor-Buneman法を用いている。
Villasenor-Buneman法は、空間メッシュの格子立方体と同じ体積の粒子で、プラズマの挙動を模擬する手法であり、粒子に一定密度の電荷を帯びさせることで、プラズマ粒子の集団特性を表し、この粒子を空間に分布させ電磁場の空間に配置することで、空間におけるプラズマ粒子の挙動が模擬できる。
また、Villasenor-Buneman法では、電流の保存法則を満たしているプラズマの場合は、空間セルの格子電流の計算によりPoisson方程式を解かずに電磁場の計算ができ、Maxwell方程式のフル解を得る事ができ、Maxwell方程式を解くには、複雑なセル境界部分を除き、超並列に配した多数の高速演算器で電流データ部分を各パイプラインに格納できるようメモリを配置することで、効率よく計算することが可能である。
また、格子電流を得るには、粒子移動後に各粒子の周囲のセルに対する、電流影響度をセル壁の貫通量によって判断するが、この判断は３次元の場合、１２以上のセル壁面を貫通する可能性があるため、それら全て内部分岐構造を再帰的手法によって処理する必要があり、汎用高速ＣＰＵのアーキテクチャである、ベクトル・キャッシュ手法を利用することはできない。
プラズマ粒子シミュレーションであるVillasenor-Buneman手法は洗練されており、且つ、電磁場のフル解が得られ、極めて優れている（非特許文献１参照）。
Jone Villasenor, "Rigorous charge conservation for local electromagnetic field solvers" Conputer Physics Communications 69(1992)306-316 North-Holland

しかしながら、非特許文献１に記載のVillasenor-Bunemanの手法においては、内部分岐構造を再帰的手法により処理する必要があるため、汎用の高速計算のベクトルキャッシュ手法が使用できない等、計算所要時間があまりにも膨大となり、現実的な時間内で計算処理できないという問題が存在する。
そこで、本発明は、プラズマ粒子シミュレーションVillasenor-Buneman法の８０％を占める計算ボトルネック部を高速に計算することができるプラズマ物理の計算機シミュレーション計算方法を提供することを目的とする。

本願請求項１に記載のプラズマ粒子シミュレーション計算方法は、空間メッシュの格子立方体と同じ体積の粒子を想定し、この粒子に一定密度の電荷を帯びさせ電磁場の空間に配置することでプラズマの挙動を模擬する、Villasenor-Buneman（ヴィラセノウ−ブーネマン）法を汎用計算機を用いたプラズマ物理の計算機シミュレーションにおいて、Newton-lorentz（ニュートン−ローレンツ）の運動方程式（式（１））の計算を、電流密度の計算（式（２））、Maxwell（マックスウェル）の方程式の計算（式（３））、粒子に働く力の計算（式（４））、と段階的に計算をするに際し、前記式（２）の計算を実行させるための計算ユニットを前記汎用計算機内に実装して、プラズマ粒子の挙動を模擬することを特徴とする。
本願請求項２に記載のプラズマ粒子シミュレーション計算方法は、前記請求項１の場合、
前記電流密度の計算をする式（２）において、初期入力データとして、ｎ個の粒子データの、粒子位置（ｘ，ｙ，ｚ）、粒子速度（ｕ，ｖ，ｗ）、電荷量（ｑ）、格子点数（ｎｘ，ｎｙ，ｎｚ）を含んだデータを与え、出力データとして、格子点上の電流成分（ｌｘ，ｌｙ，ｌｚ）を得ることを特徴とする。

本発明のプラズマ物理の計算機シミュレーション計算方法は、汎用計算機に計算専用ボードとしての計算ユニットを組み込むことにより、プラズマ物理の計算機シミュレーション計算のクリティカルセクションをハードウェアに分散化を図り、トータルとして、プラズマ物理の計算機シミュレーションの計算スピードを向上させることができる。

まず、本発明のプラズマ物理の計算機シミュレーションにおいて用いるVillasenor-Buneman法を説明する。Villasenor-Buneman法プラズマ粒子シミュレーションは、下記のVillasenor-Bunemanの基本方程式（下記の式（１）参照）を基にして、空間メッシュの格子立方体と同じ体積の粒子を想定し、この粒子に一定密度の電荷を帯びさせ電磁場の空間に配置することでプラズマの挙動を模擬する方法である。

ここで、ρ＝電荷密度を表し、Ｊ＝電流密度を表す。電流が１次元の方向に流れる場合を想定すると、Villasenor-Bunemanの基本方程式（式（１））は、電荷密度ρ、電流密度Ｊに対して、
と表せる。
ここで、電流は、一定の大きさを持った電荷素片がある面を単位時間内にどれだけ通過したのかによって定義されるので、連続の式を２精度Simplectic差分表示すると、時刻をｔ、格子座標位置をx(格子番号をi)として、式（２）は、
となる。

Maxwell方程式を同じく２次精度Simplctic差分化して電流項の時間進展は、
となる。これを、前記式（３）に代入すると、
となる。つまり、式（５）を変形すると、
となり、Poisson方程式の２次精度Simplectic差分であり、
を満たすことがわかる。

実際にこれを計算するためには、下記のように２つの場合に分けなければならない。すなわち、１次元の電荷素片が通過する面が１つの場合と、通過する面が２つの場合とである。図１は１次元の電荷素片が通過する面が１つの場合の概略説明図である。図２は１次元の電荷素片が通過する面が２つの場合の概略説明図である。
図１の場合には、
として、
により電流値を評価することができる。
図２の場合には、
として電流値を評価することができる。

次に、電流が２次元の方向に流れる場合を想定して、代表的な通過面数となる４、７、１０について、通過面数の考え方を図３乃至図５に示す。図３は２次元において電荷素片が通過する面が最小である４の場合の概略説明図である。図４は２次元において電荷素片が通過する面が７の場合、図５は２次元において電荷素片が通過する面が１０の場合である。
図３乃至図５からわかるように、２次元において電荷素片が通過する通過面数は、各象限に対して３パターン合計１２パターンの分岐構造を考えなければならない。

さらに、電流が３次元の方向に流れる場合を想定して、３次元において電荷素片が通過する面の通過面数を考えた処理方法においては、分岐構造のパターン数が多くなり急激に複雑化する。そこで、Villasenor-Buneman法では、各軸方向に垂直な面をいくつ通過するのかを判定しながら再帰処理する方法を採用する（図６参照）。ここで、再帰処理とは、電流データ部分を各パイプラインに格納できるようメモリ配置し、超並行処理する方法である。これは、電荷素片の粒子についてみれば完全に並列処理が可能であるからである。図６は、Villasenor-Buneman法における３次元再帰処理方法の概念図である。

図６から分かるように、各粒子が移動してΔｔの間に空間格子の各面を貫通する時の格子電流を求めるには、多数の分岐計算が必要になる。この結果、Villasenor-Buneman法における３次元再帰処理を、ステップを追って順に実行しようとすれば、汎用高速コンピュターの超並列演算部を用いても、複雑な分岐演算によって高速演算性能を効率よく実施することができない。

そこで、本発明においては、空間メッシュの格子立方体と同じ体積の粒子を想定し、この粒子に一定密度の電荷を帯びさせ電磁場の空間に配置することでプラズマの挙動を模擬する、Villasenor-Buneman（ヴィラセノウ−ブーネマン）法を汎用計算機を用いたプラズマ物理の計算機シミュレーションにおいて、Newton-lorentz（ニュートン−ローレンツ）の運動方程式（式（１））の計算を、電流密度の計算式（２）、Maxwell（マックスウェル）の方程式の計算式（３）、粒子に働く力の計算式（４）、と段階的に計算をするに際し、前記式（２）を実行させるための計算ユニットを前記汎用計算機内に実装して、プラズマ粒子の挙動を模擬するシミュレーション計算をすることとした。
すなわち、本発明では、格子電流の計算（電流密度の計算式（２））を、専用のハードウェア計算ユニットに委ね、他の部分の計算は汎用計算機部分で計算するシステム構成を採用し、トータルとして計算速度を向上させる構造とした。
なお、専用の計算ユニットの高速計算性能は、汎用計算機のデータ転送レートによって決定されるため、専用の計算ユニットと汎用計算機との接続は、現在ＰＣとのデータやり取りの最速仕様であるPCI-Express x16インターフェースを採用することが望ましい。

図７は、プラズマ粒子シミュレーションにおける、専用の計算ユニットの動作／実装イメージ図である。図７に示すように、汎用計算機内に、計算ユニットを実装し、電流密度の計算をする第１ステップの計算を、専用のハードウェア計算ユニットに委ね、他の部分の計算は汎用計算機部分で計算するシステム構成を採用し、有機的に融合させることで、プラズマ粒子シミュレーションにおける総合計算速度を向上させる構造とした。

図８は、空間メッシュの格子立方体と同じ体積の粒子を想定し、この粒子に一定密度の電荷を帯びさせ電磁場の空間に配置することでプラズマの挙動を模擬する、Villasenor-Buneman（ヴィラセノウ−ブーネマン）法を汎用計算機を用いたプラズマ物理の計算機シミュレーションにおいて、汎用計算機に計算ユニットを実装した場合の、Villasenor-Buneman法プラズマ粒子シミュレーション計算システム図である。
図８に示すように、まず、各粒子の運動状態は、Newton-Lorentz運動方程式に従って、式（１）→式（２）→式（３）→式（４）というようにタイムステップに順次移動する。この一連の計算処理においては、式（４）で「粒子に働く力の計算」を実行するが、式（４）中において、Ｅ（電場）、Ｂ（磁場）の計算が必要である。このためには、その前の式（３）での「マクスウェル方程式」を実行し、Ｅ、Ｂを求めなければならない。マクスウェル方程式に、Ｅ、Ｂを与えるためには、その前の式（２）で、電流保存法則によってＪ、ρを求めなければならない。

図８に示すような一連の計算を実行するにあたり、式（１）、式（３）、式（４）の各式での計算においては、粒子分の計算が殆ど独立で実施でき個別の並列計算が可能である。しかし、式（２）の電流密度の計算は、粒子が格子壁面を貫通する方向によって複数の壁に関わるケースが現れ、多数の複雑な分岐演算が発生する。この結果、汎用高速ＣＰＵのアーキテクチャである、ベクトル・キャッシュ手法等は、非効率で使用できず、シミュレーション計算の８０％の時間が式（２）の電流密度の計算に占められる状況になることが分かった。よって、本発明においては、この式（２）の格子電流計算（Villasenor-Buneman法）部分を、汎用計算機に任せず、計算ユニットに委ね、汎用計算機と有機的に融合させることで、プラズマ粒子シミュレーションの高速化を実現する。

以下に、本発明のプラズマ粒子シミュレーション計算方法において、式（２）の格子電流計算をさせるために用いたハードウェア計算ユニットの一実施例を下記に示す。ＰＣＩインターフェスとして、PCI-Express x16（転送レート約４ＧＢＰＳ）を採用した。一回の最大処理のデータ量は、５MBである。
ハードウェア計算ユニットを用いて式（２）の計算をさせるにあたって、式（１）の計算結果から出力されハードウェア計算ユニットに与えた入力データは以下のものである。
・粒子位置（ベクトル）：(x, y, z) 64bits real・・・現在位置、
・粒子速度（ベクトル）：(u, v, w) 64bits real、
・電荷量（スカラ）：ｑ 64bits real、
・格子点数：(nx, ny, nz) 32bits integer、
・出力か否かのフラッグ：
f: clear 場データのリセット、
keep 場データの保持、
back 場データの読み取り、
ハードウェア計算ユニットから計算結果として出力されるデータとして、
・格子点上の電流成分（ベクトル）：(Ix, Iy, Iz) 64bits real、
・粒子の生存状態（設定領域に存在か否か）、
を得た。この結果を、式（３）の計算をするために汎用計算機に渡し、汎用計算機においては、引き続き式（３）、式（４）の計算を継続する。
なお、ハードウェア計算ユニットにおけるデータ量／粒子は、
１個の粒子の入力データ量＝現在位置＋粒子速度＋電荷量＋格子数＋FLAG
＝ 8Bytes×3 + 8Bytes×3 + 8Bytes + 4Bytes×3 + 1Byte ＝ 69Bytes、となり、
データ転送時間（入出２回）は、5MB x 69Bytes x 2/4Gps ≒ 0.17s となる。

図９は、上記実施例の計算ユニットの内部構造を示す概略ブロック図である。

図１０は、上記実施例の計算ユニットの内部計算／処理の流れ図である。図１０に示すように、ステップ１では、前処理として入力データから計算に必要な補助データ群の生成、及び、条件分岐レジスタの設定を行う。ステップ１に要する演算時間は５００ｎｓ以下である。ステップ２では、各軸方向に通過面を１または２とする個別処理を行う。ステップ２に要する演算時間は１１００ｎｓ以下である。ステップ３では、３次元配列への補間演算・スムージング処理をを行う。ステップ３に要する演算時間は１０００ｎｓ以下である。ステップ４〜６の最終段では、補間演算／スムージング処理データとＤＤＲメモリとの通信が１５００ｎｓ以下で行われる。
総合的に、プラズマ粒子シミュレーション加速ユニットでの計算プロセスに入ってから、最終段まで１粒子あたり４１００ｎｓ以下で処理が可能となる。

図１１は、実施例の計算ユニットの計算流れ図である。図１１に示すように、多数のパイプラインで構成され２個のＦＰＧＡ内に実装されている。

図１２は、実施例の計算ユニットのハードウエアブロック図である。図１２に示すように、本実施例において、計算アーキテクチャは２個のＦＰＧＡ内に実装しており、３つ目のＦＰＧＡは追加オプションとして構成した。

図１３は、汎用計算機に実施例の計算ユニットを実装して、プラズマ粒子の挙動を模擬するプラズマ粒子シミュレーション計算方法を実行した際のシステム運用フローチャートである。

本実施例において、Villasenor-Buneman法の約８０％を占める計算ボトルネック部である、クリティカルループをハードウエア化し、大容量高速ＦＰＧＡ（Field Programmable Gate Array）内に、超並列に多数の演算器を配し、計算処理、及び、分岐処理部分まで搭載し、ハードウェアとしてボード化し、計算ユニットとして専用化したことにより、スーパーコンピュータを用いても１００時間程度必要とするシミュレーションを、約１／５程度２０時間以内に短縮することができ、優れた効果が確認できた。

本発明のプラズマ粒子シミュレーション計算方法は、今まで膨大な計算時間を要したプラズマ物理の計算機シミュレーションを、短時間でしかも研究者のディスクトップにおいて実施できる、極めて手軽で高速な専用型プラズマ粒子シミュレーション計算システムとして構築可能となり、専門的なシミュレーションの世界レベルの、最先端研究開発を支え、研究者のみならず、製造業の発展にも貢献できる。また、プラズマ物理の計算機シミュレーション以外の計算機シミュレーションにも適用可能である。

１次元の電荷素片が通過する面が１つの場合の概略説明図である。１次元の電荷素片が通過する面が２つの場合の概略説明図である。２次元において電荷素片が通過する面が最小である４の場合の概略説明図である。２次元において電荷素片が通過する面が７の場合の概略説明図である。２次元において電荷素片が通過する面が１０の場合の概略説明図である。 Villasenor-Buneman法における３次元再帰処理方法の概念図である。プラズマ粒子シミュレーションにおける、専用の計算ユニットの動作／実装イメージ図である。汎用計算機に計算ユニットを実装した場合の、Villasenor-Buneman法プラズマ粒子シミュレーション計算システム図である。実施例の計算ユニットの内部構造を示す概略ブロック図である。実施例の計算ユニットの内部計算／処理の流れ図である。実施例の計算ユニットの計算流れ図である。実施例の計算ユニットのハードウエアブロック図である。汎用計算機に実施例の計算ユニットを実装して、プラズマ粒子の挙動を模擬するプラズマ粒子シミュレーション計算方法を実行した際のシステム運用フローチャートである。

Claims

空間メッシュの格子立方体と同じ体積の粒子を想定し、この粒子に一定密度の電荷を帯びさせ電磁場の空間に配置することでプラズマの挙動を模擬する、Villasenor-Buneman（ヴィラセノウ−ブーネマン）法を汎用計算機を用いたプラズマ物理の計算機シミュレーションにおいて、
Newton-lorentz（ニュートン−ローレンツ）の運動方程式（式（１））の計算を、
電流密度の計算（式（２））、Maxwell（マックスウェル）の方程式の計算（式（３））、粒子に働く力の計算（式（４））、と段階的に計算をするに際し、
前記式（２）の計算を実行させるための計算ユニットを前記汎用計算機内に実装して、プラズマ粒子の挙動を模擬することを特徴とするプラズマ粒子シミュレーション計算方法。
前記電流密度の計算をする式（２）において、初期入力データとして、ｎ個の粒子データの、粒子位置（ｘ，ｙ，ｚ）、粒子速度（ｕ，ｖ，ｗ）、電荷量（ｑ）、格子点数（ｎｘ，ｎｙ，ｎｚ）を含んだデータを与え、出力データとして、格子点上の電流成分（ｌｘ，ｌｙ，ｌｚ）を得ることを特徴とする請求項１に記載のプラズマ粒子シミュレーション計算方法。