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

11-NA25

陰的時間積分法による核融合プラズマの非線形 MHD シミュレーションの高速化に関する研究
佐藤雅彦（核融合科学研究所）
概要	  陰的時間積分法はクーラン条件により時間刻み幅が制限されない利点があるが、連立非線形方程式を反復的に解く過程で現れる連立一次方程式を高速に処理することが求められる。本研究では完全陰解法を用いた 2 次元 MHD コード及び半陰解法による 3次元 MHD コードに対して、自動チューニング機能付き Xabclib を組み込むことにより、連立一次方程式の高速化を図った。性能評価は核融合科学研究所に設置されたHITACHI SR16000/VL1 を利用し、SMT 実行を含む 128 スレッド実行環境において行った。本研究の性能評価試験によれば、反復解法ライブラリ Lis を組み込んだ MHD コードに対して、Xabclib を組み込んだ MHD コードは 10 倍〜100 倍と十分に高速であることが確認され、自動チューニング機能の導入による高速化を確認できた。

1. 研究の目的と意義	 れるため、より高解像度の格子が必要とされるよう	  核融合炉を実現するためには、高温高密度のプラになる。高解像度格子を使用すると、陽的な時間積ズマを安定に閉じ込める必要があり、プラズマ中で分法ではさらに小さい時間刻み幅が要求される。こ発生する磁気流体力学(MHD)不安定性の非線形現象のため、シミュレーションで使用されているプラズの解明が重要である。このため、核融合科学研究所

マ抵抗値は、実際の実験パラメータよりも 3 桁も大では、一流体 MHD モデルに基づいた MHD 非線形きい値を使用しているのが現状となっている。

シミュレーションコード(MIPS コード)を開発し、	  そこで本研究では、陰解法を用いた MHD シミュ

MHD 不安定性の非線形現象の解明に向けた研究をレーションコードの開発することにより、これらの進めている。MIPS コードは 4 次のルンゲクッタ法課題を解決することを試みる。陰的時間積分法を使による陽的な時間積分法を使用しているため、時間用すれば、クーラン条件に制限されることがないた

積分における時間の刻み幅がクーラン条件により制め、高解像度の格子を使用し、かつ、時間の刻み幅限されるという問題点がある。MHD 不安定性の中にを大きくしたシミュレーションが可能となる。これ

おいて興味のある現象は、プラズマ抵抗値に依存すにより、実験結果の解釈や予測を可能とする、実験る抵抗性 MHD 不安定性であり、その非線形現象のに即したより現実的なパラメータでの核融合 MHD

時間スケールはミリ秒のオーダーである。これに対シミュレーションが可能となり、核融合研究の進展し、プラズマコア領域においては、一流体 MHD モに大きく貢献できると考えられる。

デルに含まれる最も短い時間スケールの現象は、マ	  陰的時間積分法においては、連立非線形方程式をイクロ秒のオーダーを持つ速い磁気音波である。こニュートン法により解く必要がある。ニュートン法のため、現在使用している陽的な時間積分法では、では反復計算により解を求めることになるが、その時間刻み幅が速い磁気音波によって制限されること反復過程において大規模な連立一次方程式を解かなになる。特に、抵抗性 MHD 不安定性の成長率はプければならない。このため、陰的解法の高速化にはラズマ抵抗値に依存し、プラズマ抵抗値が小さくな大規模疎行列に対する高速な線形ソルバーの開発がるとともに成長率も小さくなる性質があるため、プ必須となる。本研究では、自動チューニング機能付

ラズマ抵抗値が小さい計算に対しては計算時間を長き数値計算ライブラリ Xabclib を導入することによくとらなければならなくなる。また、プラズマ抵抗り、大規模連立一次方程式の高速な処理を目指す。

が小さくなるとより小さな空間スケールの構造が現

この大規模疎行列問題は前処理付きクリロフ法によ1

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

り解くことが有効であると考えられる。一般的な問れてきた。実行時にアルゴリズム（実装法）選択

題と同様に、クリロフ法では有効な前処理法の導入のため、入力データとノード内スレッド数の変化が必要であるが、プラズマが持つ強い非線形性のたに応じた AT 機能を提供することを目的としている。

めに大規模疎行列の条件数が非常に大きいことが問OpenMP で並列化しているので、ノード内高スレッ題である。そのため、MHD 非線形シミュレーション

ド実行時の AT 効果が、Xabclib では重要となる。	 に適した前処理方法の研究を行い、大規模並列計算	 機の性能を有効に活用した高速な線形ソルバーの適	  本研究課題は昨年 H22 年度から開始されたもので用を行っていくことも本研究の課題である。ニューある。H22 年度においては、完全陰解法コードの基トン法で用いられる高速な線形ソルバーは、他の核礎的検証のために 2 次元モデルに対する完全陰解法融合プラズマシミュレーションコードにおいても有コードの開発を行った。この 2 次元コードにおいて用である。固有値問題による MHD 線形安定性コー

は、Lis を用いて連立一次方程式を並列処理している。

ド、他のスキーム（有限要素法、コンパクトスキー

このコードを用いて抵抗拡散時間が依存するプラズム等）による MHD 非線形シミュレーションコード

マ抵抗値、及び、アルフヴェン波の伝播速度が依存等へも適用することが可能であろう。特に、MHD 線するプラズマ密度のパラメータと完全陰解法コード

形安定性コードがプラズマの実時間制御に使用されの収束までの反復回数との関係を調べ、完全陰解法る場合、その計算時間の高速化が重要な課題となっコードの有効性を確認した。これと平行して、Xablibている。したがって、本研究の成果はより幅広いプに実装された並列性とキャッシュ再利用性を高めラズマシミュレーション研究の発展に貢献できるとた Blocked	 BSS	 (Branchless	 Segmented	 Scan)	 法

期待される。の性能検証を最大で 1 ノードあたり 128 スレッド	  Xabclib は、東京大学情報基盤センターの片桐を実行が可能な HITACHI	 SR16000/VL1（以降、SR16K）中心としたグループにより開発が進められているを用いて行った。これにより、高スレッド実行で自動チューニング（以降、AT）機能付きライブラ

Xabclib の AT 機能の有効性を確認できた。

リである。数値計算ライブラリにおける AT 技術は	  本年度においては、昨年度の研究成果を基に２密行列や FFT など信号解析ライブラリで成功を収次元完全陰解法コード及び 3 次元半陰解法コード

めてきた。これは、階層キャッシュに代表されるに Xabclib を組み込み、Xabclib による性能向上ハードウェアの複雑化と数値計算ライブラリ特有の評価を行った。の複雑なチューニング手法による開発コストの増	 大が背景にある。	 	 

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

一方，プラズマシミュレーションを含む多くの数値シミュレーションでは、疎行列に対する連立(1) 共同研究を実施した大学名と研究体制	 一次方程式の求解が行われる。そのため、疎行列

東京大学	 用の数値計算ライブラリ（疎行列ライブラリ）の(2) 共同研究分野	 性能が重要となる。疎行列ライブラリは、入力デ超大規模数値計算系応用分野	 ータの特徴を考慮した実装をしないと、多くの場(3) 当公募型共同研究ならではという事項など	 合、高性能が達成できない。そのため、入力デー高並列スレッド環境として、核融合研究所が	 	 	 タの特徴を自動抽出し、最適実装を自動選択する

所有する HITACHI	 SR16000/VL1 を利用し、	 方式が研究されている．	 

SMT 実行を含む 128 スレッド実行環境での	 このような状況の中、我々は疎行列ライブラリにXabclib を組み込んだ陰解法 MHD コードの性	 おける実行時 AT の機能開発を目的として AT 機能能評価を行った。	 

付き数値計算ライブラリ Xabclib の開発が進めら	 2

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

3. 研究成果の詳細と当初計画の達成状況	 !f1=(#fi +2 + 8fi +1 # 8fi#1 + fi#2 )!x12"x	 !2 f
1=(#fi +2 + 16fi#1 # 30fi + 16fi#1 # fi#2 )!x 212"x 2(1) 研究成果の詳細について	 	  本研究の最終的な目標は 3 次元完全陰解法コードの開発であるが、本年度においては、昨年度に開発した 2 次元モデルでの完全陰解法コード及び、 陽解法版 MIPS コードでは、4 次精度の陽的ルンゲ半陰解法を用いた 3 次元コードに対して Xabclib

クッタ法を用いて時間積分を行っている。本コーを導入して性能評価を行った。本節では開発のベドは MPI により並列化され、核融合科学研究所設ースとなる MIPS コードの概要、3 次元半陰解法コ置のプラズマシミュレータにて、MHD 不安定性のードの概要、2 次元完全陰解法コードの概要につ非線形挙動の解明に使用されている。	 

いて述べた後、Xabclib による性能評価について	 述べる。	 

Z	 !

3.1	 MIPS コ ー ド の 概 要 	 
MIPS は、下記の一流体モデルの MHD 方程式を円柱座標系 (R, !, z) のもとで解くコードである。	 

!!
= "#·(!v)!t% v 2 (*!v!
= "!" $ v " !# ''' ** " #p + j$ B 	 '& 2 *)!tR

4+ #[#!(#·v)]" # $(#!")
3!B= "# $ E!t!p= "#·(pv) "(! "1)p#·v!t4+ (! "1)["#$ 2 + "#(#·v)23	 + % j·(j " jeq )]	 E = "v$ B + %(j " jeq )法コードの開発を行ってきた。ファラデー則とオj=	 
図 1.	 MIPS コードで使用されている円柱座標系	 

3.2	 3 次 元 半 陰 解 法 コ ー ド の 概 要 	 	  MIPS コードの陰解法化の第一段階として、オーム則における拡散項のみを陰的に取り扱う半陰解

1#$Bµ0ーム則より、磁場の時間発展方程式は、	 !B= " #(v# B) $ "! # j + " #(! jeq ) + !"2B!t$ = #$v
ここで、! はプラズマ密度、v はプラズマの速度、とかけ、 ! を拡散係数とした磁場の拡散項が現れ

p はプラズマの圧力、 B は磁場、 E は電場、 j はる。この式を分離し、第一ステップにおいて	 プラズマ電流密度を表している。jeq は平衡での電!B= " #(v# B) $ "! # j + " #(! jeq ) 	 !t

流密度である。また、! はプラズマの抵抗値、! はを陽的に時間積分を行った後、	 粘性、 ! は比熱比である。4 次の中心有限差分により空間方向の離散化を行い、1 階微分、2 階微分!B= !"2B 	 !tはそれぞれ、次式の様に離散化されている。	 

を陰解法により時間積分を行う。	 
3

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究 最終報告書 2012 年 5 月	  実際のプラズマにおいては、プラズマ中心から運動方程式における粘性項及び圧力方程式におけ周辺に向かうにつれプラズマの温度は減少するたる散逸からのエネルギー取得項については無視しめ、プラズマ抵抗が周辺に従うにつれて大きくないている。	 

る。このため、プラズマコア領域ではアルフヴェ	  本研究の性能評価試験においては、初期平衡とン波の伝播時間によりクーラン条件が決まるが、して磁場は y 成分のみを持つとし、その分布を	 周辺部においては磁場の抵抗拡散時間でクーラン条件が決まるようになる。上記の半陰解法コード

By = B0 tanh(x / d) 	 

により抵抗拡散時間によるクーラン条件の制限をとして与えた。ここで、 d は磁場のシアーの幅を緩和でき、より実際に近い抵抗分布のもとでシミ表すパラメータである。圧力分布は、平衡方程式ュレーションを行えるようになる。	 

!p = j" B から決定している。今考えている問	 

題では、	 

3.3	 完 全 陰 解 法 を 用 い た 2 次 元 MHD コ ー ド の概 要 	 dBdp= ! y By 	 dxdx	  完全陰解法コードの基礎的検証のため、図 2 にから求められる。また、初期速度は 0 と仮定する。

示すような 2 次元スラブ配位を対象とした 2 次元

この初期平衡のもとで磁場に微少な摂動成分を与MHD コードを開発している。(x,y,z)座標系のもとえ、物理変数の時間発展を計算していくことになで、物理変数は z 方向に一様と仮定する。(x,y)

る。	 

y x ：固定境界座標系において、y 方向に周期境界条件、x 方向にラブモデルの対応は、トーラスプラズマにおけるy ：周期境界z ：一様半径方向とスラブモデルの x 方向が対応し、トーBは固定境界条件を考える。トーラスプラズマとスラスプラズマにおけるトロイダル角(φ)方向がス

x

ラブモデルにおける z 方向に対応していることになる。スラブモデルにおいて、簡単のために速度、磁場の z 方向の成分は 0 とすれば、(x,y)座標のもとで MHD 方程式は、	 

# !2!(!vx ) !(!vy )!!
! 2 &(=""+ D %%% 2 + 2 (( ! + S !
%$ !x!t!x!y!y ('!vx!t!vy!t!Bx!t!By図 2.	 2 次元スラブモデル	 	 = "(vx!
!
1 !P 1+ vy)vx "" JB!x!y! !x ! z y	  この 2 次元コードでは、時間積分法として完全= "(vx!
!
1 !P 1+ vy)vy "+ J z Bx!x!y! !y!

連立非線形方程式を解く必要があり、ニュートン=陰解法を組み込んでいる。陰的時間積分法では、法により解を反復的に求める必要がある。連立非!
!
(vx By " vy Bx ) "["(J z "Jeq )]!y!y線形方程式を	 F(u) = 0 	 とおくと、推定解の近傍で線形化した方程式	 !
!
(vx By " vy Bx ) +["(J z "Jeq )]!t!x!x!Fdu = "F(u) 	 
22 &#!(Pv)!(Pv)!P!
! ((y%x!u=""+ # %% 2 + 2 ( P + SP!t!x!y!y ('$% !x	  から解の更新量 du をもとめ、	 ="	 u ! u + du 	 と書き下すことができる。ただし簡単のために、により推定解を更新する。この過程を F(u) が要求
4

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

精度を満たすまで反復する。記憶容量の削減等を

図 3 に 2D 問題の実行時間[秒]を載せる。図１で目的として、ヤコビアン行列を明示的に与えないは、スレッド数を 16、32、64、128 と変化させて方法が主に行われているが、本研究で開発した 2いる。	 

次元 MHD コードでは、ヤコビアン行列は解析的に与えている。	 	 
3.3	 Xabclib に よ る ソ ル バ 高 速 化 	 
3.3.1	 性 能 評 価 試 験 	 	  本節では、Xabclib2011 利用による疎行列反復解法ソルバ（GMRES(m)法を用いたもの）を、2D および 3D の陰解法 MHD コードに適用した場合のシミュレーション実行時間を評価する。評価環境を以	 下にまとめる。	 
l

計算機環境：核融合研に設置された	 

図 3	 2D 問題の実行時間[秒]。N=153,984。

HITACHI	 SR16000/VL1 の 1 ノード（64 コア）。	 
l

利用スレッド数：16、32、64、128(SMT)	 

l

評価対象	 
Ø

NNZ=3,648,064。要求誤差：1.0d-12。	 
前処理は「無し」を指定。	 
図 3 に よ り 、 Lis に 対 し Xabclib2011( 以代表者佐藤による陰解法 MHD コード。2D

降、”Xabclib”と記載)は約 24 倍高速化されてい問題、および 3D 問題のプログラムで評価る。収束までの反復回数は、Lis は 10 回〜11 回、する。	 

Xabclib では 6 回であった。	 

²

利用言語：Fortran90	 	 ²

コンパイラ：日立最適化 Fortran90

3.3.3	 3D 問 題 の 評 価 結 果 	 

コンパイラ。コンパイラオプション	  図 4 に、3D 問題の実行時間[秒]を載せる。	 

は”-Oss	 ‒omp”。	 
Ø

比較対象のライブラリの詳細	 
²

Xabclib2011 ： OpenATLib2011のOpenATI̲LINERSOLVE 関数を呼びだす。
GMRES(m)法を利用。数値計算ポリシは「実行時間」。疎行列‐ベクトル積の実装、リスタート周期(m)の調整には、自動チューニングを適用。	 
²

²

Lis	  ver.1.2.26 に お け る 関 数
lis̲solve() を 呼 び 出 す 。 GMRES(m)	 

法を利用。リスタート周期(m)は、デ

図 4	  3D 問題の実行時間[秒]。N=245,760。

フォルトである 40 を指定。	 

NNZ=3,356,160。要求誤差：1.0d-6。前処理はLis,	 Xabclib2011 とも、ライブラリ「ILU(0)」を指定。32 スレッドを利用。	 

や行列の初期化などの実装は同等の	  図 4 では、X 軸にシミュレーションにおけるスユーザ関数として実現。	 

テップ(kstep と記載)を取っており、kstep ごとに	 反復解法ソルバが呼ばれる。	 

3.3.2	 2D 問 題 の 評 価 結 果 	 

図 4 から、Xabclib ではシミュレーションが進5

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

むにつれ、実行時間が減少している。しかし Lis
では、実行時間は変化しない。Xabclib を利用することで Lis に対し、シミュレーションの開始時に約 12 倍の高速化、シミュレーションが進むにつれて約 98 倍の高速化が達成される。	 	  この 3D シミュレーションでは、各 kstep において前回の解を利用することで、ソルバ内の収束加速を行う実装になっている。したがって、理論的にはシミュレーションが進むにつれて実行時間は減少するはずである。Lis および Xabclib では、	 前回の解ベクトルを利用するように、実行時オプ

図 6	  3D 問題における 1 反復あたりの実行時間

ションの指定、および初期ベクトル設定（Xabclib）[秒]。前処理時間を含んだ時間。N=245,760。

を行っている。したがって、Lis の実行時間は予NNZ=3,356,160。要求誤差：1.0d-6。前処理は想外である。そこで、収束するまでの反復回数を「ILU(0)」を指定。32 スレッドを利用。	 

調べた。図 5 にそれを載せる。	 	  図 6 より、kstep=20 以降では Xabclib は反復 0回なため、ほとんど前処理時間になると考えられる。この場合には、Lis より Xabclib は 1 反復あたりの実行時間が約 3.7 倍高速化されている。一方、kstep=20 未満では、Lis より最大でも約 1.9倍高速化されるに留まる。	 	 
3.3.4	 考 察 	 	  2D 問題、3D 問題ともに、Lis に対して Xabclibは 10 倍〜100 倍と十分に高速であることは疑いよ	 
図 5	  3D 問題における収束するまでの	 うがない。したがって、自動チューニング機能の反復回数。N=245,760。NNZ=3,356,160。	 

導入により高速化が実現されたといえる。	 

要求誤差：1.0d-6。前処理は「ILU(0)」を指定。	  2D 問題の結果を考察する。前処理なしの GMRES32 スレッドを利用。	 

実行については、Lis は 10 回〜11 回反復であり、	  図 5 より、Lis では一定回数の 25 回反復をしてXabclib では 6 回反復と少ない。これは、GMRES 法いるのに対し、Xabclib では 4 回、2 回、1 回と、の実装の違いにより生じた可能性が高い。具体的kstep が進むたび減少している。kstep=20 以降でには、Xabclib では直交化部分に「完全 QR 分解」

は反復は 0 回（前処理のあと１回だけ反復解法のを使っている。このことで演算の高効率化と、場処理を行い収束したということ）となる。したが合により高精度化されることが予想される。結果って、Lis の本バージョンにおいては、解ベクトとしてそれが、高速化と反復回数の原書に貢献しルの再利用に関しするバグの可能性が高い。	 

た可能性がある。	 	  1 回あたりの反復回数を調べることで、より詳	  3D 問題の結果を考察する。前処理のみの実行時細な実行時間評価をすることができる。そこで、

で Xabclib は Lis に対してより高速である。この

図 6 に 1 反復あたりの実行時間[秒]を載せる。た

理由は、ILU(0)の実装に関して、Xabclib と Lis

だし反復回数 0 回は 1 回とみなしている。	 の実装の違いにより生じたと考えている。Xabclibの ILU(0)実装では、対角部分のみしか値更新をし6

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究 最終報告書 2012 年 5 月#1%%%% !
%% "t%
D1= %%% 0%%%%%%% 0%$

ない軽量化した ILU(0)実装を採用している。このことで、容易に収束する行列では、Lis より ILU(0)前処理が高速化される可能性が高い。	 	  以上は、実装の違いによる定説的な考察といえる。問題の性質やスレッド数を変化させた状況下ことが必要である。これは今後の課題である。	 	 

01!
"t0

#%%1 )%% !
%% !0 )x
L=%%%%% ! 1 )%% !0 )y$

で、性能モニタを用いて定量的に性能解析をする&(((#((1%%((%% !
(("t0 ((, D2 = %%((%%(%%$ 0
1 (((!
("t (('0B0 )!0 )y
J0!0(&((((((1 (((!
("t ('0#&())%%(%% ! (P0 *) ! (P0 *) ((B0 ) (&((()y%% )x((! !
((%!0 !0 )x ((((((, U = %%% ) (B *)0((0%% )y(((((((%%0(((()%%(('!
(B*)
0((%%$ )x 0('J0	  昨年度実施した陰解法を用いた 2 次元 MHD コー# n+1 (&# n (&%% P%% P ((# n+1 &(# n &%% n+1 ((%% vx (1 %% n (((1 %% vx (((((%% Bx ((, r2 = !
u = %% Bx ((, v = %% n+1 (, r1 = !
%("t %% n (("t %%% vyn (((%% n+1 ((%% vy ((($'$'((%%$ By ('%%$ By ('ドの性能評価では、プラズマ密度が減少するにつとおけば、	 れて反復法の収束が悪化し、計算時間が大きく増! D U $&!
!
$## 1&&# u $&& = ## r1 &&& 	 ##&&###&&& ###" L D2 &%" v % #" r2 &&%3.4	 physics-based	 preconditioning の導入	 

大することがわかった。これは、密度が減少するにつれてアルフヴェン波の伝播速度が増加し、双	 曲型問題の性質が強くなったためと考えられる。と書くことができる。このとき、Schur 分解によこの問題を解決するために、物理変数に対してり下記のように分解できることが知られている。	 

Schur 分解を行うことで双曲型問題が放物型問題! D## 1###" Lに変換し、収束を向上させる手法が提案されている。この手法は MHD 問題においても有効であることが確認されている（L.Chacon	 et	 al,	 J.	 Comput.	 

ここで、	 

Phys.,	 178,	 15	 (2002)）。本研究ではその手法のS = D2 - LD1-1U 	 効果を確認するために、下記の 2 次元線形 MHD モデルに対して評価を行った。	 
!vx!t!vy!t!Bx!t#By="="=である。これを利用すると、	 

!Bx &((1 !P1
1 #% !B" J 0By ! %% y "(!0 !x !0!0 %%$ !x!t (('1 !P1+ JB!0 !y !0 0 x

U $&& !# I0 $&&!# D1 0 $&&!## I D1-1U $&&&& = ##&##&&&#-1
D2 &%& ##" LD1 I &&%#" 0 S &&%##" 0I &%	 (1) u = D r(2) v = r - Lu 	 (3 ) v = S v(4) u = u - D Uv
*-1
1 1*

*

2-1 *	 !
(B v ) + !!"Bx!y 0 x*-1
1#= $ (B0vx ) + !!"By
#t
#x
#P
#
"= $ (P0vx ) !
(P v ) + !!"P#t#x
"y 0 yから解を求めることができる。いま、	 ここで、下付の 0 は平衡量を表している。本テスであるから、	 

D1!1 = !"t I 	 トではプラズマ密度の値と収束性能に注目するため、非拡散項のみを陰的に処理し、その後で拡散

S = D2 + !t LU 	 項を陽的に処理する。このとき、陰的に処理をすとなる。この行列 S はもとの行列と比べて次元がる非拡散項の部分は、	 小さくなるとともに対角優位となるため、クリロフ法の収束が早くなることが期待できる。計算コードの性能評価においては、図 2 で与えられている平衡に対して行った。プラズマ密度分布は	 
ρ0=(1-ρa)(1-x2)4+ρa と仮定した。周辺での密度7

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

値ρa と、時間ステップを 20 ステップ進めるのに(a)グローバルな行列問題を解かない手法	 

要した実計算時間の関係を図 7 に示す。計算では単純に連立一次方程式をプロセスごとに分割した

16 コアを使用し、Lis ライブラリに組み込まれて際、となりのプロセスが担当する変数の値をどのいる GMRES 法を用いた。従来の Schur 分解を使用ように取り扱うかが問題となる。もっとも簡単なしないコードと比較して、改良コードでは計算時方法は、隣のプロセスが担当する変数は固定境界間が大幅に削減できていることが確認できる。特としてそのプロセスで取り扱う方法である。このに、周辺密度が小さくなるにつれてその違いが大ようにすることで、各プロセスで独立に連立一次きく現れており、physics-based	 preconditioning方程式を解くことができる。連立方程式を解いたの導入が効果的であることを確認できた。	 

後、各プロセスの境界部分の値を更新し、再び連立一次方程式を解く。この作業を境界部分の値が実計算時間（s）収束するまで反復させる。この手法ではグローバ
ルな行列問題を解く必要がない、プログラミンが行いやすいという利点があるものの、連立一次方程式を反復手法により解く上で、収束性能が悪化する問題点がある。	 
(b)静的縮約法	 
Schur 分解を用いた領域分割を行うことより計算時間の削減が期待できる。しかしながら、MPI 並列を用いた Schur 行列問題の高速解法、および、周辺密度ρa	 	  疎行列同士の積の高速化などの課題がある。この

図 7. 周辺密度に対する計算時間の関係

2 手法のうち、まず(a)の手法により MPI+OpenMP(2)	 当初計画の達成状況について	 ハイブリッド並列コードの開発を進め、既存ライ最終的な目標は 3 次元非線形コードの陰解法化でブラリと Xabclib との性能比較試験を進めていく。

あったが、現在のところ 2 次元非線形コードの陰その次の段階として、より複雑な(b)の手法を用い解法化まで達成している。2 次元モデルのもとでた計算コードを今後開発していく予定である。	 

ソルバーの改良及び検証を進め、その有効性を確

5. 研究成果リスト	 認した。ここで行った成果をふまえて 3 次元コー(1)	 学術論文（投稿中のものは「投稿中」と明記）	 

ドへの拡張を進めていく予定である。	 

(2)	 国際会議プロシーディングス	 

4. 今後の展望	 (3)	 国際会議発表	 	  これまでに開発してきた陰解法を用いた MHD コ(4)	 国内会議発表	 

ードは、OpenMP を使用したノード内並列を行って	  ・佐藤雅彦、プラズマシミュレータシンポジウ	  	 いる。さらなる高速化のため、また、高解像度メ	  	  ム 2011、自然科学研究機構	  核融合科学研究	 

ッシュを用いた大規模計算を行うためには MPI を	  	  所主催、2011	 年	 9	 月	 13 日(火)~14	 日	 (水)	 

用いたノード間並列も行っていく必要がある。陰	  「陰的時間積分法を用いた MHD	 コード開発」	 

解法コードにおいては、前処理付きクリロフ法に・片桐孝洋、T2K シンポジウム 2011＠京大、	 より連立一次方程式を解く必要があるため、大規	 	 2011 年 11 月 30 日（水）	 

模ノード間並列にて計算速度の向上を図るために「e-Science プロジェクト「シームレス高生産・	 は領域分割法の導入が必要不可欠と考えられる。	 高性能プログラミング環境」成果報告」	 

領域分割法を用いた連立一次方程式の解法として、 (5)	 その他（特許，プレス発表，著書等）	 
下記の 2	 つの方法について検討を行う予定である。	 
8