J. Plasma Fusion Res. Vol.90, No.9 (2014)５６２‐５６７

研究論文

実験データと理論計算データを融合した統合輸送解析システムの開発
Development of the Integrated Transport Analysis System
Unifying Experiment and Simulation Data
江 本 雅 彦，鈴 木 千 尋，鈴 木 康 浩，横 山 雅 之，關

良 輔，居 田 克 巳

EMOTO Masahiko, SUZUKI Chihiro, SUZUKI Yasuhiro, YOKOYAMA Masayuki,
SEKI Ryosuke, and IDA Katsumi
自然科学研究機構 核融合科学研究所（原稿受付：２０１
４年２月３日／原稿受理：２０１
４年７月７日）

実験データと理論計算データを融合した統合輸送解析システムを開発した．このシステムにより，従来，別物として取り扱われてきた，LHD の実験データと理論計算データを統合化し，ユーザは両者を区別することなくデータを取り扱うことができるようになった．このシステムにより，平衡計算および，その平衡に基づく NBI加熱分布評価，エネルギー・運動量バランス解析に至る一連の解析作業をパッケージ化することが可能となり，
LHD 実験適用型統合輸送解析コード TASK3D の運用に大きな貢献を果たし，原理的に LHD の全放電，電子温度分布計測がなされた全てのタイミングに対してエネルギー・運動量バランス解析を実施することができるようになった．
Keywords:
transport analysis, unified system, numerical simulation, database, TSMAP

格化トロイダルフラックス %を用いて "!%で表され，##$$

１．序論

は実効プラズマ小半径で，プラズマ体積 " とプラズマ幾何

核融合プラズマを記述する統合シミュレーションコード

"中心 !!を用いると，"!"の関係にある．
!!!!##$$の開発は，国内の TASK3D［１］や TOPICS-IB［２］だけでな
く，EU の JETTO［３］や米国の TRANSP［４］等，世界中で行

一方，LHD はプラズマ形状が複雑であり，測定されるプ

われており，予測シミュレーションだけでなく実験解析にラズマ断面は測定器毎に異なるため（図１），規格化座標も用いられている．LHD 実験に対応したシミュレーション
コードの開発では，シミュレーション側からの視点に重点が置かれており，これらのコードに精通していない研究者が，実験データの解析を行うのは不便であった．例えば，
輸送解析ではプラズマパラメタが一定であるとみなせる時間領域のみを対象とした静的な解析が主として行われ，プ
ラズマ実験で得られた全プラズマ放電時間を対象とした動的な解析を行うには，必要な全ての実験データをシミュレーションコードが必要とするデータ形式に変換する必要があり，このような解析はあまり行われてこなかった．
この理由の一つは，実験データとシミュレーションコードが必要とするデータ形式は異なっており，シミュレーションコードを実行するには，必要とされる複数の収集データから必要な形に変換する必要があるため，非常に手間がかかるためである．例えば，TASK3D では，入力パラ

図１メタとして，径方向を実効半径 $#$$で表した座標系（以下，規格化座標系）で与える必要がある．ここで実効半径 $#$$"で定義される．また，"は規格化小半径で，規は$#$$!##$$

計測毎に異なるプラズマ断面：YAGレーザを用いたトムソン散乱計測では，プラズマの横長断面を計測しているのに対し，CH3OH 遠赤外レーザを使った計測ではプラズマの縦長断面を測定している．

National Institute for Fusion Science, Toki, GIFU 509-5292, Japan
Corresponding author’s e-mail: emoto.masahiko@nifs.ac.jp

５６２!2014 The Japan Society of Plasma
Science and Nuclear Fusion ResearchContributed Paper

Development of the Integrated Transport Analysis System Unifying Experiment and Simulation Data

M. Emoto et al.

系への変換作業が必要となるが（図２），これには，プラズのはデジタイザから取得された生のデータのため，計測しマの平衡を求め，全計測で統一的で一貫性のある変換作業た研究者以外には，物理的に意味のあるデータに変換するが必要である．

ことは困難である．LHD で収集された全ての物理量のデー

また，シミュレーション実行には複数の物理データを必

タを統一的なインタフェイスを用いて参照できるようにす要とするが，物理データ間には依存関係があり，ある物理るという目的で開発されたのが，解析サーバシステムであデータが更新された場合には，依存する物理データの更新

る［７］．このシステムでは，スプレッドシート等，多様なが必要となる．後述する TSMAP［５］では平衡を決定するツールを利用して読み込むことができるよう，物理量に変ため，トムソン散乱計測による電子温度分布を利用してい換したデータを解析データフォーマットというカンマ区切るが，例えば，暫定値として登録されていたトムソン散乱りのテキストファイルとして提供している．解析サーバシ計測による温度分布の値が，その後の校正により更新されステムは，図３に示すようにファイルサーバとリレーショた場合は，平衡が変わるため，規格化座標系での様々な物

ナル・データベース（RDB）から構成されており，ファイ理分布のデータの更新が必要となり，数値計算で利用するルの保存場所，ファイルサイズ，登録時刻等のメタ情報を全てのデータを更新する必要がある．

RDB で管理し，ファイルそのものはファイルサーバに保存

本研究は実験データと数値計算を統合化し，実験データしている．また，クライアント・サーバ間の通信は gsoap

からシミュレーションコードが必要とするデータを自動作［８］を 利 用 し た，HTTP に よ る SOAP
［９］を 使 用 し て い成し，全実験データでの数値計算を行えるような枠組みをる．このシステムを利用することで，LHD 実験で収集また作るものである．このためには，１．物理データアクセスへ

は解析を行った様々な結果を，解析サーバシステムのクラの改善，２．座標変換，３．依存データの自動更新が必要とイアントプログラムやライブラリを利用することで統一さなる．次章では，これらの機能を実現するために開発したれた方法で取得できるようになった．

システムについて説明を行う．

２．
２ 座標変換の必要性
解析サーバを使うことで，物理データの統一化を行うこ

２．実験データと数値計算の統合とができるようになったが，これらは個々の計測器が観測

２．
１ 物理データの提供している物理座標における物理量であり，輸送解析等の数

LHD 実験では１００を超え る 計 測 装 置 が 取 り 付 けられ，値計算で利用するには，規格化座標系での物理量に変換す

LABCOM システム［６］によってデータ収集および提供がる必要がある．また，Thomson 散乱計測では，電子密度の行われている．ところが，このシステムに登録されている分布を得ることが出来るが，これを絶対値に変換するにはFIR 計測による電子密度の線積分の値によって絶対値校正する必要があり，このために座標変換が必要となる．これは，LHD ではプラズマ形状が複雑であり，観測しているプラズマ断面は計測毎に異なるためである．例えば，Thomson 散乱計測ではプラズマの横長断面を測定しているのに対し，FIR 干渉計ではプラズマの縦長断面を測定している．
２．
２．
１ 平衡データベース
実座標から規格化座標系への変換を行うために，VMEC［１０］コードを用いて平衡データベースの作成を行った．
VMEC コードは，プラズマ中の全エネルギーを変分原理に

図２

図３

座標変換：（a）
長半径方向の測定と（b）
短半径方向の測定から絶対値校正を行うには（c）のように規格化座標系へのマッピングを行う必要がある．

５６３

解析サーバ概要：解析サーバはメタデータを格納する RDBとファイルサー バ か ら な る．ク ラ イ ア ン ト と の 通 信 はSOAP を利用して行う．

Journal of Plasma and Fusion Research Vol.90, No.9 September 2014

ことができた．

より最小化することで，MHD 平衡を求めるものである．
このデータベースは，ある磁場配位とプラズマ形状におけこのデータベースを利用するために，PV-Wave,FOR-る磁気面形状を磁気座標系 $%"&"$%で表し，そのフーリエTRAN,C 等の各種コンピュータ言語から呼び出せるライブ係数を格納している．但し，%は規格化されたトロイダルラリを提供している．ユーザはこのライブラリを利用し，フラックス（プラズマ境界で %#"），&はトロイダル方向，上記７つのパラメタと実座標の配列を渡すと，実座標に対$はポロイダル方向を表す．

応する規格化座標の配列を得ることができる．検索処理のヘリカルコイルが作る磁場配位は３つの物理量，Rax：

実装は先の解析サーバシステムをベースとした，クライア

磁気軸位置，Bq：ヘリカルコイルとポロイダルコイルの四ント・サーバ方式を採用しており，検索処理はすべてサー重極成分の比，%：ヘリカルコイルのピッチパラメタ，で決バ側で行っている．

定できる．ただし，%#!#"で，!，$はコイルの周$"#*#期数と極数，#*は小半径，"は大半径である．また，VMEC存しており，一方，入力パラメタ等のメタ情報は RDB に保の計算にはプラズマ圧力および電流の分布をトロイダルフ存している．サーバはクライアントから受け取ったパラメ

ラックス %の多項式として与える必要があるが，計算の高タに一致する入力パラメタをRDBから検索し，一致するも速化のため，幾つかの分布について前もって磁気面形状をのがみつかれば，ファイルの内容をクライアントに返す．

計算しておき，実際の実験で得られた分布を用いる場合完全に一致するパラメタが存在しない場合は，パラメタ空は，あらかじめ計算した近傍の値から補間して求めるよう間における近傍の値を用いて，補間したデータを返す．必にした．ここで，プラズマ圧力については表１にあげる５

要とする座標の多くは既存の計測器が観測を行っている視種類の分布関数を使用した．また，プラズマ電流について線方向の座標なので，変換の高速化のため，主要な計測器は，ホローな電流分布も考慮し，表２に示す１１種類を使用の視線方向の座標については，あらかじめ計算しておき，VMEC による計算結果はサーバ側にファイルとして保した．また，それぞれの分布関数の形状をピーキングファ

その結果を視線データベースとして利用している．

クターという１つのパラメタで表すことで，探索を高速に２．
２．
２ 最適な平衡の決定行っている．ただし，ここで用いたピーキングファクター

前節で述べたデータベースを利用して，プラズマにおけは s の全領域で積分し平均した値と中心値の比である．

る 平 衡 を 求 め る プ ロ グ ラ ム TSMAP［６］を 作 成 し た．

全ての磁場配位の組み合わせで変換テーブルを作るに

TSMAP はプラズマのピークベータ値と圧力ピーキング

は，表３に示す計７つのパラメタ空間上で行う必要があフ ァ ク タ ー を 変 更 し な が ら，与 え ら れ た 磁 場 配 位 で，り，１つの磁場配位あたり，７８６５（＝１１×５×１３×１１）通Thomson 散乱計測で得られた電子温度分布を規格化座標りの計算が必要となる．データベースの作成にあたり，実系へマッピングし，プラズマの内側の電子温度と外側の温験でよく使う磁場配位から計算を始めており，２０１３年末現度との差が最小になる形状パラメタを探索し，そのプラズ在で２５パターン，約２０万点の計算が終了している．計算量

マ に お け る 平 衡 を 求 め る プ ロ グ ラ ム で あ る（図４）．

が 膨大 な た め，２００９年 当 時，大 型 汎 用 機 を 用 い た 計 算Thomson 散乱のデータを利用するのは，この計測が常時が，１磁場配位あたり，２０〜３０日かかっていた．（ただし，

表２

これは計算速度が遅いためではなく，むしろ，ユーザが多電流分布関数と電流ピーキングファクター（lpf）．
電流分布関数

くジョブの実行時間の割り当てが少ないことによる．
）

$$%
"!%
#
$%
"!%
$%"!%

この計算を高速化するために，PC６台，計２８CPU コアを用いて並列に計算させることにより，５日に計算を短縮する#

表１$"!%%%$%$'$"!%"!% %!!!
$%#%$)$"!%"!% %!!!
$%#%
$
#$
"!%
"!% %!"!
$%#%$&$"!%"!% %!"!
$%#%$($"!%"!% %!"!
$%#
%
$
"$"!%
"!% %!#!
$%#
%
$%$
"!%
"!% %!#!
%#

圧力分布関数と圧力ピーキングファクター（pf）．
圧力分布関数
$
"!% %
$%
$
"!%
"!%%%

pf
1.41
2.14

$%
"!%
$
$%
"!%
%
$%
"!%

3.00
4.00
5.00

%#

#

表３パラメタ
Rax
Bq%
p0
pf
IpIpf

平衡データベースパラメタ一覧．

磁気軸位置
四重極磁場
ピッチパラメタ説明範囲
3.50−4.00−50〜200
1.129〜1.262

分割数
１１６８

単位
m％
−

ピークベータ値
圧力ピーキングファクター
1 T あたりのトロイダル電流値
電流ピーキングファクター

0〜10
1.41〜5.00−150〜150−12.60〜4.00

１１５１３１
１％−kA/T−５６４Ipf
4.0
3.0
2.0
1.41
0.71
0.21
−0.49
−1.49
−3.06
−5.91
−12.60Contributed Paper

図４

Development of the Integrated Transport Analysis System Unifying Experiment and Simulation Data

M. Emoto et al.

TSMAP による平衡の決定：規格化座標系にマッピングしたプラズマの内側と外側の電子温度分布の差が最小になるよう，パラメタを変えながら探索を行う．

行われており，かつ空間計測点が多くマッピングの精度が期待できるからである．ただし，トロイダル電流値はロゴスキーコイルの値から決定され，また，最終的な結果に対する電流ピーキングファクターの影響は小さいため，2.0に固定している．これによって，プラズマの平衡が決まり，実座標から規格化座標系へのマッピングができるようになった．ま た，こ れに よ り FIR の 測 定 デ ー タ を 使 っ てThomson 散乱で得られた電子密度の絶対値校正を行うことが可能となった．
しかしながら，典型的なショットでは，Thomson 散乱計
測は１ショットあたり，８００フレーム以上の計測を行って図５おり，ライブラリの合計呼び出し回数は１０，０００回に達する．そのため，平衡計算には１０分以上かかる．一方，LHDの短パルス実験では３分周期で放電が行われるため，平衡並列計算によるクライアント側の高速化：時間フレーム毎にプロセスを生成し，複数の CPU を使って並列計算を行うことで高速化を行う．

計算が実験サイクルに追いつかない．この問題に対処する様に記述することができるが，並列処理を書くことのできため，当初６台のクライアント PC を用意し，各 PC が１る言語で，このような記述が可能な言語ということから，ショットずつ計算を行っていた．これによって，１ショッ

Python［１１］を選択した．Python では，numpy という配列

トの計算速度そのものは速くならないが，数ショット遅れ処理パッケージを利用することで，PV-Wave に似た記述で，実験の放電ペースに追いつくことが可能となった．しができるため，移植が比較的容易であると判断したためでかしながら，ショットの解析を一刻も速く行い，次回以降ある．また，並列化には multiprocessing というパッケージの放電パラメタの設定に反映する必要があることから，一を用いた．このパッケージを使い，時刻毎のループで行っプログラムでの高速化が必要となった．ている処理をそれぞれのプロセスに分割し，同時に実行す２．
２．
３ 平衡計算の高速化ることで計算時間の短縮を行った．

TSMAP の高速化を行うため，クライアント・サーバのこの方法により，クライアント側の並列化によりクライ双方で高速化を行った．

アント側の処理を高速化できたが，複数のクライアントかクライアント側では並列実行を行うことで高速化した．

らの同時リクエストが増大することになり，サーバ側の拡TSMAP によるによる平衡の計算では，個々の時刻におけ張も必要となった．サーバ側は，元々マルチプロセスで動る計算は独立しているため，それぞれの時刻における計算いているため，複数クライアントからの同時リクエストをは並列して行うことが可能である（図５）．したがって，マ処理することは可能であるが，同時に実行されるプロセス

ルチコアCPUの計算機を利用することで，各フレーム計算数がCPUコア数を超えることがあるため，複数のサーバをを同時に実行すれば，一放電あたりの計算時間が短縮する用意し，同時実行できるプロセス数の増加を測った．これことが可能である．フレーム毎に並列して計算を行うにには LVS（Linux Virtual Server）を利用することで仮想的は，TSMAP のプログラムをマルチスレッドあるいはマルに複数のPCでのサービスを１台のPCで行っているかのよ

チプロセスを利用し，並列計算するよう書き換える必要がうにしている（図６）．現在，仮想サーバは LVS サーバ１ある．しかしながら，オリジナルの TSMAP プログラムは台と３台の１２コア（２４スレッド同時実行可能）の実サーバ

PV-Wave で記述されており，並列プログラムを書くのにから構成されている．クライアントから LVS サーバにリクは適していない．そこで，プログラムを別の言語に移植すエストがあると，LVS サーバは最も接続数の少ない実サーることにした．PV-Wave では配列の演算を普通の変数のバにリクエストを転送し，そのサーバから直接クライアン
５６５

Journal of Plasma and Fusion Research Vol.90, No.9 September 2014

図６

LVS によるサーバ側の高速化：LVS を利用し，３台の物理PC を１台に見せかけることで，クライアント側の変更なしに増大したクライアントリクエストに対応する．

トに計算結果を返すという方式（ダイレクト・ルーティン
グ）をとることで，クライアントからのリクエスト数の増加に対応した．

図７

物理データの依存関係：輸送解析の物理データ dytrans6を 作 成 す る に は，cxmap6，TSMAP 等 が 必 要 で あ り，TSMAP を実行するには，thomson，ip，wp 等のデータが必要となる．

図８

AutoAna 概要：解析サーバに新規にデータが登録されると，ジョブキュー登録デーモンはジョブをキューに登録する．登録されたジョブは複数の実行サーバにより並列実行される．

以上のような高速化手法を用いることにより，TSMAPの実行時間は１０分あまりから２５秒にまで短縮することが可能となり，実験シーケンス中での処理が可能となった．
２．
３ 自動解析
TSMAP で行った平衡計算の結果を利用して，規格化座標系における物理量を計算するには，TSMAP による平衡データを解析サーバに登録する必要がある．放電サイクル中に迅速に解析を行うためには，TSMAP の平衡計算が終わり，平衡データが登録されると，すぐにこのデータを利用する解析プログラムを実行する必要がある．また，ある実験データが解析データとして登録された後に，データの校正やプログラムの改良等により，更新が行われることがあるが，この場合には依存するデータの更新も必要となることがある．例えば，TSMAP が実行時に利用するデータが更新された場合，TSMAP の平衡データの更新が必要となり，さらに，この平衡データに依存するデータの更新も必要となる．図７は計測データから，輸送解析に必要なデータを作成し，結果を格納するまでのデータの流れを図示したものであるが，左端の thomson のデータが更新されたとすると，AutoAna サーバは TSMAP が実行可能かどうた際には，TSMAP の再計算により平衡が変わるため，そかチェックする．TSMAP の実行には thomson，ip，wp，…

れに依存する１０以上のデータが更新される必要がある．そ…のデータが必要であり，もしこれらのデータが揃っていのため，筆者らはデータ間の依存関係を管理し，自動的にれば，TSMAP の実行ジョブを投入する．次に TSMAPデータを登録・更新する AutoAna システムの構築を行っの計算結果が解析サーバに登録されると，このデータに依た［１２］．

存 す る ermap6，cxsmap6…実 行 可 能 で あ る か ど う か処理の流れは図８に示すようになる．

チェックを行い，実行可能であれば順に解析データ登録プ

AutoAna システムは，１．
ジョブキューサーバ，２．ジョログラムの実行ジョブを登録していく．ジョブは FIFO

ブ登録デーモン，３．ジョブ実行デーモンの３つ の コン（先入れ先出し）の待ち行列で管理されており，ジョブ実行ポーネントから構成される．

デーモンは待ち行列からジョブを登録された順番に取り出

解析サーバにデータが登録されると，IP マルチキャストし，処理して行く．AutoAna サーバが提 供 し て い る のを利用してネットワーク上の PC に新たなデータが登録さは，解析プログラムを起動する枠組みだけで，解析プログれたことを通知する．ジョブ登録デーモンはこの通知を受

ラム自体は，それぞれの解析に精通した研究者が作成してけると，そのデータから作成されるべきデータがあるかどいる．これらの解析プログラムは，最低限のインタフェイうかチェックし，そのデータの登録プログラムが実行可能

スを実装していれば，研究者の必要に応じて PV-Wave やであれば，ジョブキューサーバに登録プログラムのジョブFORTRAN 等，様々な言語で作成することができる．

を投入する．

ジョブ実行デーモンは，複数のPCで動いており，同時実例えば，図７で thomson が新たに解析サーバに登録され行可能である．２０１３年度の実験中は５台のPCを使い，
４０プ
５６６

Contributed Paper

Development of the Integrated Transport Analysis System Unifying Experiment and Simulation Data

M. Emoto et al.

ロセスを同時に動かしていたが，この数は動的に変更可能である．例えば，解析プログラムの変更が行われた場合には，登録された過去のデータを全て計算し直す必要があり，場合によっては，数万〜数十万件のデータの再計算が必要となることがある．このような場合でも，一時的にCPU の空いている PC を利用し，ジョブ実行デーモンを動かすことで，急な要求にも柔軟に対応できる．

３．数値計算への貢献
本研究により，統合輸送解析ソフト TASK３D の運用に大きな貢献を果たした．例えば，図９は放電中における，密度で規格化した熱流束とイオン温度勾配の関係を２０ミリ秒毎に図示したものである［１３］．このようなパワーバランスの解析が容易に行えるようになったのは，本研究で構築したシステムにより，TSMAP に基づく平衡計算および，その平衡に基づく NBI 加熱分布評価等，実験データからエネルギー・運動量バランス解析に至る一連の作業を自動化することで，全ショットにおける上記のデータが解析サーバという単一システムよって参照が可能となり，実験データと統合シミュレーションの連携が強化されたためである．

まとめ本システムにより，従来，別物として取り扱われてきた，
LHD の実験データと数値計算を統合化した．これより，ユーザは数値計算に必要なデータを容易に入手できるよう

図９ （a）
内部輸送障壁（ITB）の外側（規格化小半径＝0.84）の領域での電子密度で規格化された熱流束とイオン温度勾配の関係 （b）
ITB の内側（規格化小半径＝0.57）での両者の関係（参考文献［１４］より引用）．になり，実験シーケンス中における実験データを入力パラメタとした，高次の輸送解析を行うことのできる枠組みを提供することが可能となった．また，TSMAP の高速化で，時間フレーム毎に分割した並列処理を行うことによる高速［３］G. Cenacchi and A. Taroni, Rapporto ENEA RT/TIB 88/5 (1988).
［４］http://w3.pppl.gov/transp/［５］C. Suzuki et al., Plasma Phys. Control. Fusion 55, 014016(2013).
［６］H. Nakanishi et al., Fusion Sci. Technol. 58, 445 (2010).
［７］M. Emoto et al., Fusion Eng. Des. 81 2019 (2006).
［８］http://www.cs.fsu.edu/〜engelen/soap.html［９］http://www.w3.org/TR/2000/NOTE-SOAP-20000508/［１０］S.P. Hirshman and O. Betancourt, J. Comput. Phys. 96, 99(1991).
［１１］http://www.python.org/［１２］M. Emoto et al., Fusion Eng. Des. 89, 758 (2014).
［１３］K. Ida et al., Phys. Rev. Lett. 111, 055001 (2013).

化を実証することができた．この並列処理を，本システム内部で使用している，他のサブシステムに適用していくことで，実験サイクル中における解析が可能であることを示せた．

謝辞

このシステムの開発の機会を与えていただき，また，貴重なアドバイスをいただいた LHD 実験の関係者の皆様方に感謝いたします．この研究は LHD 計画予算（ULHH０１４）によって行われました．

参考文献［１］M. Yokoyama et al., Plasma Fusion Res. 7, 2403011 (2012).
［２］N. Hayashi et al., Plasma Fusion Res. 6, 2403065 (2011).

５６７