GPGPU によるメッシュ法を用いた磁束クリープ
フローのパラメータ推定における計算の高速化

電子情報工学科 小田部研究室
学籍番号 10232057
鶴田祐基

平成 29 年 2 月 16 日
電子情報工学科目次
目次 ................................................................................................................................... i
序論 .................................................................................................................................. 1
1.1

はじめに........................................................................................................................ 11.2

磁束ピンニング............................................................................................................ 11.3

磁束クリープ・フローモデル .................................................................................... 21.3.1

磁束クリープによる磁界 .................................................................................... 2

1.3.2

磁束フローによる電界 ........................................................................................ 5

1.3.3

ピンポテンシャル ................................................................................................ 6

1.3.4

磁束クリープ・フローモデル ............................................................................ 8

1.4

不可逆磁界.................................................................................................................... 9

1.5

メッシュ法.................................................................................................................. 10

1.6

シンプソン法.............................................................................................................. 10

1.7

GPGPU (General-purpose computing on graphics processing units) .......................... 12

1.8

GPU の構造................................................................................................................. 13

1.8.1

ストリーミングマルチプロセッサ（SM） .................................................... 13

1.8.2

レジスタ .............................................................................................................. 14

1.8.3

共有メモリ .......................................................................................................... 14

1.8.4

テキスチャメモリ .............................................................................................. 14

1.8.5

コンスタントメモリ .......................................................................................... 14

1.8.6

グローバルメモリ .............................................................................................. 14

1.8.7

キャッシュ .......................................................................................................... 14

1.9

CUDA .......................................................................................................................... 15

1.10

GPU における並列化の階層..................................................................................... 15

1.11

実験の目的 .............................................................................................................. 15

第2章

計算手法 ....................................................................................................... 162.1

評価値𝑑の計算方法 ................................................................................................... 16

2.2

プログラムの流れ...................................................................................................... 17

2.3

実行環境...................................................................................................................... 18

2.4

解析に用いるデータ.................................................................................................. 18

第3章

結果・考察 ................................................................................................... 19
-i-

3.1

CPU と GPU の計算時間の比較 ............................................................................... 19

3.2

GPU で一度に計算する条件について ..................................................................... 22

第4章

まとめ ........................................................................................................... 25

謝辞 ................................................................................................................................ 26
参考文献 ........................................................................................................................ 27

図目次
図 1.1 磁束クリープにおける磁束バンドルとエネルギーの関係 ............................ 3
図 1.2

磁束クリープにおける磁束バンドルとエネルギーの関係 ........................ 5

図 1.3

磁束線と平行方向の相関距離と超伝導層の厚さによって変わる磁束バン

ドルの体積 ................................................................................................................... 7
図 1.4 不可逆領域と可逆領域 ........................................................................................ 9
図 1.5 2 つのパラメータを格子状に組み合わせた概略図 ....................................... 10
図 1.6 積分区間を 3 等分した場合のシンプソン法 .................................................. 11
図 1.7

GeForce GTX 960 ................................................................................................ 12

図 1.8 GPU の構造 ......................................................................................................... 13
図 2.1 距離dの模式図.................................................................................................... 16
図 2.2 プログラムの流れ .............................................................................................. 17
図 3.1 パラメータの分割数と CPU、GPU の計算時間 ............................................ 19
図 3.2 パラメータの分割数と評価値の関係 .............................................................. 20
図 3.3 分割数を変化させた場合の磁束クリープ・フローモデルによる理論値のEJ特性と実験値のE-J特性との比較 ........................................................................... 21
図 3.4

GPU で一度に計算する条件を変更した結果 ................................................. 23

表目次
表 2.1 パラメータの範囲 .............................................................................................. 16
表.2.2

実行環境 ............................................................................................................... 18

表.2.3

使用する線材の諸元 ........................................................................................... 18

表 3.1 分割数が 8、11、14 の時のdの値 .................................................................... 21
表 3.2 それぞれのプログラムにおける計算条件 ...................................................... 22- ii -

序論
1.1

はじめに
1908 年にオランダの Leiden 大学の Kamerlingh Onnes が世界で始めてヘリウムの液体化に成功した。また、彼は 1911 年に液体ヘリウムを用いて水銀の電気抵抗を測定する過程において、4.2 K 以下の条件下で電気抵抗がほぼ 0 となることを発見した。
さらに 1933 年には Fritz Walther Meissner によって超伝導体の完全反磁性（Meissner 効果）が発見されている。この 2 つの現象を超伝導現象と呼ぶ。超伝導体は電気抵抗がないため、電流を損失なく通電することができると期待されていた。しかし、超伝導体はある磁界を越えると超伝導状態を維持できないことが分かった。超伝導状態から常伝導状態に替わる温度を臨界温度𝑇c 、磁界を臨界磁界𝐵𝑐 と呼ぶ。超伝導状態になる
条件はあるものの、電気抵抗が 0 となることは有用性が高いため、現在でも研究は続けられている。
その後も超伝導現象の発現についてはしばらく不透明なままであったが、1977 年に John Bardeen，Leon Neil Cooper，John Robert Schrieffer の 3 人によって提唱された
BCS 理論によって超伝導現象のメカニズムが明らかにされた。BCS 理論によると、
𝑇c は 30 K を超えないだろうと予想されていた。しかし、1986 年に Johannes Georg
Bednorz と Karl Alex Müllar によって臨界温度が 30 K を超える酸化物超伝導体が発見された。その 1 年後には常圧下の液体窒素の沸点（77.3 K）を超える酸化物超伝導体
として、YBa2 Cu3 O7−δやBi2 Cr2 Ca2 Cu3 O10+δ (𝛿は微小量)などの酸化物超伝導体が発見された。これらは銅酸化物超伝導体もしくは高温超電導体と呼ばれ、安価な液体窒素
中で超伝導状態を示すことから応用が期待されている。
また、2001 年に青山学院大学の秋光純教授のグループがMgB2 を発見し、2008 年には東京工業大学の細野秀雄教授らが鉄ヒ素系超伝導体の発見されている。このように、現在に至るまで超伝導体の発見や研究が盛んにおこなわれている。

1.2

磁束ピンニング

超伝導体は磁気的な振る舞いによって第一種超伝導体と第二種超伝導体に分類される。第一種超伝導体の場合。超伝導状態においては常に完全反磁性を示すが、𝐵c を越えると超伝導状態が消失する。第二種超伝導体では、下部臨界磁場𝐵c1 までは第一種超伝導体と同じように完全反磁性を示すが、これ以上の磁場を加えると、超伝導状態と常伝-1-

導状態が混在する混合状態となり、上部臨界磁場𝐵c2 を超えると完全に常伝導状態になる。
第一種超伝導体の𝐵c と第二種超伝導体の𝐵c1 は𝐵c2 と比較すると非常に小さい。したがって、工業用の用途としては、𝐵c1 と𝐵c2 の間の混合状態で用いられる。混合状態では、常伝導状態の部分に磁束線が進入する。この状態で電流を流すと、磁束線は Lorentz 力
𝑭L を受ける。Lorentz 力𝑭L は超伝導体に流れる電流密度𝐽と進入した磁束密度𝐵を用いて𝑭L = 𝑱 × 𝑩と表すことができる。𝐹L を受けた磁束線が速度𝑣で動いたとすると、電磁誘導
によって、𝑬＝𝑩 × 𝒗の電界が生じる。しかし、実際には磁束線の運動を妨げるように
Lorentz 力とは反対向きに同じ大きさの力が働く。この力をピン力といって、ピン力が働く作用のことを磁束ピンニングという。また、ピンニングに作用する、常伝導析出物
や結晶粒界などの超伝導状態にならない部分のことをピンニングセンターと呼ぶ。単位
面積あたりのピン力を𝐹p とよぶ。𝐽𝐵 = |𝑱 × 𝑩|とすると、𝐽𝐵が𝐹p の最大値を超えなければ電界は発生しない。したがって、電気抵抗 0 で流すことのできる最大の電流密度𝐽c は、
𝐽c ＝𝐹p /𝐵と表すことができる。𝐽c を臨界電流密度といい、超伝導体の応用において、𝑇cや𝐵c2 と同じくらい重要なパラメータである。𝐹p を向上させることにより、𝐽c は増加する。ピンニングセンターは超伝導体の不均一さ、つまり欠陥であり、人工的に不純物を加えてもピンニングセンターとして機能する。ピンニングセンターを人工的に導入することで、より大きな𝐽c を得ることが知られ、人工ピン導入についての研究がなされてい
る[1]。

1.3

磁束クリープ・フローモデル

1.3.1 磁束クリープによる磁界
磁束クリープとは、ピンニングセンターによって固定されている磁束線が熱振動によって確率的にピンニングセンターから外れてしまう磁束線の運動のことである。理論的には外部環境が変化しない限り、超伝導電流は減衰しないと考えられるが、実際に超伝導電流を長時間測定すると減衰している。外部環境が変化していないにもかかわらず、
超伝導電流が時間によって減衰しており、ピンニングによる超伝導電流が恒常的なものでないことを示している。また、高温になると熱運動が盛んになるため超伝導電流の減衰が著しくなり、高温超伝導体の場合には𝐽c が 0 になることもある。
超伝導体には磁束線が量子化して進入するため、磁束クリープの際にも、磁束線は束になって移動すると考えられている。この磁束線の束を磁束バンドルという。電流が流れている状態における一つの磁束バンドルについて考える。この時の磁束バンドルの変化は図 1.1 のように表すことができる。磁束線は Lorentz 力を受けているため、図のように右下がりの図のようになる。磁束線がピン止めされた状態から抜け出すためには、-2-

エネルギーバリアの𝑈や𝑈′を超えなければならない。熱振動がなければ、図の状態で安
定となる。

𝑈’
𝑈

flux Bundle

図 0.1

𝑈0

𝑎f

磁束クリープにおける磁束バンドルとエネルギーの関係

磁束バンドルは温度𝑇の下では、熱エネルギー𝑘B 𝑇（𝑘B はボルツマン定数）に依存した確率で𝑈を超える。この確率は、Arrhenius の式より、exp(−𝑈/𝑘B 𝑇)で与えられる。超伝導体に進入した磁束線は、格子間隔af の三角格子を作るので、磁束バンドルはaf 移動
することで次の安定状態に移動すると考えられる。したがって、磁束バンドルの熱振動数を𝑣0 とおくと、この時の Lorentz 力方向の平均の磁束線の移動速度𝑣+ は、

𝑣+ = 𝑎f ν0 exp (−

𝑈
)
𝑘B 𝑇

(1.1)

となる。Lorentz 力の反対方向の磁束バンドルの平均の移動速度も考慮すると、全体での移動速度𝑣は、

𝑣 = 𝑎f ν0 [exp (−𝑈𝑈′) − exp (−
)]
𝑘B 𝑇
𝑘B 𝑇

(1.2)

となる。𝑣0 とaf はそれぞれ、

𝜈0 =𝜁𝜌f 𝐽c02π𝑎f 𝐵

(1.3)

1

2𝜑0 2
𝑎f = (
)
√3𝐵

(1.4)

で与えられる。ただし、𝜁は、ピンニングセンターの形状によって決まる定数であり、点状のピンニングセンターの場合はζ ≅ 2π、ピンニングセンターがaf 以上の非超伝導粒子の場合は、ζ = 4となることが知られている。𝜌f はフロー抵抗であり、𝐽c0 は磁束クリープがないと仮定した場合の仮想的な臨界電流密度である。
この𝐽c0 は、経験的に求められており、

-3-

𝑚

𝑇 TCN_T𝐵 TCN_B
𝐽c0 = 𝐴 (1 − ( )) 𝐵 𝛾−1 (1 −
)
𝑇𝑐
𝐵c2

(1.5)

と表現されている。ここで、𝐴と𝛾はピンニングパラメータと呼ばれるパラメータ、
TCN_T と TCN_B はスケール則内での次数である。例えば、銅酸化物超伝導体においては、TCN_T＝2、TCN_B＝2 である。したがって、𝑬＝𝑩 × 𝒗の関係と(1,1)式より、生じきる電界の大きさは、𝐸 = 𝐵𝑎f ν0 [exp (−𝑈𝑈′) − exp (−
)]
𝑘B 𝑇
𝑘B 𝑇

(1.6)

となる。すなわち、磁束クリープによって式(1,6)で与えられる電界が生じていることになる。
磁束バンドルの中心位置𝑥に対するエネルギー変化は、一般的には図 1.1 であらわされる正弦洗濯板状ポテンシャルになり、

𝐹(𝑥) =𝑈0sin(𝑘𝑥) − 𝑓𝑥
2

(1.7)

で近似的に与えられる。ただし、𝑈0 /2はポテンシャルの振幅、𝑘 = 2𝜋/𝑎f は波数、𝑓 = 𝐽𝐵𝑉は Lorentz 力である。磁束バンドルが平衡状態にある時を𝑥 = −𝑥0 とすると、𝑥 = 𝑥0 でエネルギーは極大となる。したがって、それぞれの位置でのエネルギーの変化量は 0 となる（𝐹′(𝑥) =0）から、𝑥0 は、

𝑥0 =

𝑎f
𝑓𝑎f
cos −1 (
)
2π
𝑈0 π

(1.8)

と求まる。また、𝑈 = 𝐹(𝑥0 ) − 𝐹(−𝑥0 )であるから、𝑈は、

𝑓𝑎f
𝑓𝑎f
𝑓𝑎f
𝑈 = 𝑈0 sin[cos−1 () ]−cos−1 (
)
𝑈0 π
π
𝑈0 π
1
2 2

2𝑓
2𝑓
2𝑓
= 𝑈0 [{1 − () } −cos −1 (
)]
𝑈0 𝑘
𝑈0 𝑘
𝑈0 𝑘

(1,9)

と表される。ただし、sin⁡(cos−1 𝑥) = √1 − 𝑥 2 を用いた。熱振動がないと仮定すると、𝑈 =0の理想的な臨界状態となると考えられる。そのための条件は、

2𝑓
2𝐽c0 𝐵𝑉
=
=1
𝑈0 𝑘
𝑈0 𝑘
を満たすことである。したがって、

2𝑓
𝐽
=
=𝑗
𝑈0 𝑘 𝐽c0(1.10)

したがって、(1,9)式は(1,10)式を使うと、
1

𝑈(𝑗) = 𝑈0 [(1 − 𝑗 2 )2 − 𝑗 cos−1 (𝑗)]-4-

(1.11)

また、𝑘 = 2π/𝑎f と、(1,10)式より、

𝑈 ′ (𝑗) ≅ 𝑈 + 𝑓𝑎f = 𝑈 + π𝑈0 𝑗(1.12)となる。(1.12)式を利用して、(1.6)式を整理すると、

𝐸 = 𝐵𝑎f 𝜈0 exp (−𝑈
π𝑈0 𝑗) [1 − exp (−
)]
𝑘B 𝑇
𝑘B 𝑇

(1.13)となり、磁束クリープによる電界を求めることができる。

1.3.2 磁束フローによる電界
磁束クリープの状態からさらに電流を流すと、ピン力が Lorentz 力を支えきれなくな
り、磁束線が連続的に運動している状態になる。これを磁束フローと呼ぶ。磁束フローの状態のとき、図 1.2 に示すようにポテンシャルは変異に対して単調減少する。

𝑈’

図 0.2

𝑈0

磁束クリープにおける磁束バンドルとエネルギーの関係

磁束フローで発生する電界のみを評価するため、磁束クリープがない状態を考える。
超伝導体に電流が流れ、外部磁界が印加されている時、単位面積当たりの磁束線に働く
Lorentz 力は𝑱 × 𝑩で与えられる。また、1.3.1 節で述べたように、磁束線は Lorentz 力と逆方向にピンニング力を受ける。Lorentz 力の方向の単位ベクトルを𝛅 = 𝒗/|𝒗|とすると、
つりあいの式は、

𝑱 × 𝑩 − 𝜹𝐹p = 0(1.14)となる。この時の電流密度は、𝐽c0となることから、|𝑱| = 𝐹p0 /𝐵 = 𝐽C0 の関係が得られる。
磁束フローの影響を考えると、粘性力が働くため、-5-

𝑱 × 𝑩 − 𝜹𝐹p −𝐵
𝜂𝑣 = 0𝜙0(1.15)となる。ただし、ηは粘性係数である。𝐽c0 ＝𝐹p0 /𝐵と、𝐸＝𝐵 × 𝑣を使って、𝐽について解
くと、

𝐽＝𝐽c0 +

𝐸
𝜌f

(1.16)となる。ここで、𝜌f = 𝐵𝜙0 /𝜂はフロー比抵抗である。また、𝜌f は、常伝導状態の抵抗密
度𝜌0 を用いて、𝜌f = 𝜌0 𝑇𝐵/𝑇c 𝐵c2で求まる。(1.16)式を𝐸について解くと、

𝐸 = 𝜌f (𝐽 − 𝐽c0 )(1.17)となり、磁束フローで発生する電界が求まる。

1.3.3 ピンポテンシャル
ピンニングセンターの持つピンポテンシャルは、磁束クリープによる超伝導電流の緩和や不可逆曲線を決定する際に重要となる。この、ピンポテンシャルのエネルギー𝑈0 は、
̂0 と、磁束バンドルの体積𝑉の積で表す
単位体積あたりの平均化したピンポテンシャル𝑈ことができ、

̂0 𝑉
𝑈0 = 𝑈(1.18)

̂0 は、Labusch パラメータαL と、相互作用距離𝑑i を用いて、となる。この𝑈

̂0 =
𝑈

1
𝛼 𝑑2
2 L i

(1.19)と表すことができる。ピン力密度𝐹p と𝐽c0の関係は、

𝐹p = 𝐽c0 𝐵⁡ ⁡ ⁡ ⁡ ⁡ = αL 𝑑i(1.20)

また、𝑑i は、𝑑i = af /𝜁の関係があるので(1.20)式とこれらの関係から、𝑈0 は、

𝑈0 =1
𝐽 𝐵𝑎 𝑉2𝜁 c0 f

(1.21)と表すことができる。この式から、𝑈0 は、ピン力の強さだけでなく、磁束バンドルの体積にも左右されることが分かる。
磁束バンドルはある領域内で秩序が保たれた磁束線の集団であり、運動する際には、集団で移動する。この秩序は磁束格子のピンニング相互作用によってもたらされており、磁束バンドルの大きさはピンニング相関距離によって決定される。磁束線と平行方向および垂直方向の相関距離をそれぞれ、𝐿、𝑅とするとこの 2 つの相関距離によって、磁束バンドルの体積が求められる。また、磁束バンドルの体積が求まると、式（1.21）から、理論的な𝑈0 の値を求めることができる。この磁束バンドルの体積を模式的に表した図を
図 1.3 に示す。

-6-

𝑑

𝐿

𝑑

𝐿
𝑩

𝑅
𝑅
𝑅

𝐿<𝑑

図 0.3

𝑅
𝐿<𝑑

磁束線と平行方向の相関距離と超伝導層の厚さによって変わる磁束バンドルの体積

𝑅は磁束格子間距離af と等しいかその実数倍と考えられており、(1.22)

𝑅 = 𝑔𝑎f22であらわす。𝑔 は磁束バンドル中の磁束線数である。𝑔 は、熱力学的な原理から磁束クリープ下での𝐽c が最大となるように決定されると考えられている。その式は、

𝑔2 = 𝑔e2 [5𝑘B 𝑇𝐵𝑎f ν0
log (
)]
2𝑈e
𝐸c

(1.23)

であらわされる。ただし、𝑔e2 は、完全格子を組んだ際の𝑔2 の値であり、𝑈e は、𝑔2 = 𝑔e2のときの𝑈0 の値である。また、𝑔e2 は、

𝑔e2 =

0
𝐶66
2π𝐽c0 𝐵𝑎f

(1.24)

0で与えられる。𝐶66
は完全格子の剪断定数であり、

0
𝐶66
=

𝐵c2 𝐵
𝐵
(1 −
)
4𝜇𝐵c2
𝐵c2

(1.25)

であらわされる。ただし、𝐵c は熱力学的磁場である。
一方で、𝐿は弾性理論より、
1

1

𝐶44 2
𝐵𝑎f 2
𝐿=( ) =()
𝛼L
𝜁𝜇0 𝐽c0(1.26)となる。ここで、𝐶44は磁束線の曲げ歪みに対する弾性定数であり、𝐶44 = B 2 /𝜇0 で与えられる。𝐿は、超伝導体が薄い場合にはその厚さ𝑑に制限がかかる場合があり、バンドル
サイズは、𝐿 < 𝑑と、𝐿 < 𝑑の場合でそれぞれ異なる。𝐿 < 𝑑の場合、ピンニングは三次元-7-ピンニングとなり、図 1.3 のように平行方向のバンドルサイズは𝐿となる。したがって、
𝑈0 は、
1/2

𝑈0 =

0.835𝑔2 𝑘B 𝐽c0
𝜁 3/2 𝐵1/4

(1.27)となる。次に、𝐿 < 𝑑の場合には、2 次元ピンニングとなり平行方向のバンドルサイズ
は、𝑑となる。よって、𝑈0 は、

𝑈0 =

4.23𝑔2 𝑘B 𝐽c0 𝑑
𝜁𝐵1/2(1.28)

ただし、7/4(2/√3)(𝜙07 /𝜇02 )1/4
≈ 0.835𝑘B
2
3

3

2 2
( ) 𝜙02√3≈ 4.23𝑘B
2の数値的関係を用いた。ここで、𝑘B を用いているのは、𝑈0 を熱エネルギー𝑘B 𝑇と関連つけるためである。

1.3.4 磁束クリープ・フローモデル
以上に述べてきたように、超伝導体で発生する電界については、磁束クリープによるものと磁束フローによるものの 2 つがある。この 2 つを考慮して超伝導体に発生する電界を理論付けする計算モデルが磁束クリープ・フローモデルである。1.3.1 項、1.3.2 項から、磁束クリープによる電界𝐸cr 、磁束フローによる電界𝐸ff は、
𝑗 ≤ 1の範囲では、

𝑈(𝑗)
π𝑈0 𝑗
) [1 − exp (−
)]
𝑘B 𝑇
𝑘B T

(1.29)

𝐸ff = 0(1.30)

𝐸cr = 𝐵𝑎f ν0 exp (−𝑗 > 1の範囲では、

𝐸cr = 𝐵𝑎f ν0 [1 − exp (−

π𝑈0
)]
𝑘B T

𝐸ff = 𝜌f (𝐽 − 𝐽c0 )(1.31)(1.32)と表すことができる。ただし、(1.32)式は𝐸cr が残留するためにこのような式であると仮定されている。これらの 2 つによって生じる電界𝐸′を
2 + 𝐸2𝐸 ′ = √𝐸cr
ff

(1.33)とする。
酸化物超伝導体では、超伝導体内部は不均一であるために、ピン力は分布していると-8-考えられる。したがって、ピン力の分布は正規分布に従うと考えて、(log 𝐴 − log 𝐴m )2]
2𝜎 2

𝑓(𝐴) = 𝐾exp [−(1.34)

で表現することにする。ここで、𝐾は規格化定数、𝜎 2 はピン力の分散、𝐴m はピン力の最頻値を表す。この𝐴の分布を考慮すると、発生する電界は、∞𝐸(𝐽) = ∫ 𝐸 ′ 𝑓(𝐴)d𝐴(1.35)

0となり、この式を計算することで、理論的な𝐸-⁡ 𝐽特性が得られる[2]。

不可逆磁界

1.4

第二種超伝導体は、第一種超伝導体と比較して、超伝導状態を維持することができる。
第二種超伝導体は、混合状態では超伝導部分が残っているために上部臨界磁場𝐵𝑐2 までに抵抗 0 で流すことができる電流が存在すると考えられていた。実際には、𝐵c2 以下であってもピンニングが働かなくなり𝐽c が 0 となる。この𝐽c = 0と𝐽c ≠ 0の境界の磁界のことを不可逆磁界𝐵i と呼び、温度ごとの不可逆磁界の値をつないだ曲線のことを不可逆曲線と呼ぶ。これらをまとめて、図 1.4 に示す。

𝐵
𝐵c2 (𝑇)

𝐵i (𝑇)

可逆領域
𝐽c = 0不可逆領域
𝐽c ≠ 0𝑇
図 0.4

不可逆領域と可逆領域

-9-

1.5

メッシュ法

メッシュ法とは、パラメータの範囲を等間隔で分割し、分割した複数のパラメータを組み合わせることで多次元の格子状にし、パラメータの全ての組み合わせを計算する手法である。計算の精度を上げる場合には、メッシュの目を細かくすればいいため、分割数を大きくするだけで実装できる。ただし、分割数を大きくすると、計算時間も増大す
る欠点がある。具体的には、𝑛種類のパラメータ全てを𝑘分割すると、(𝑘 + 1)𝑛 回の計算が必要となる[3]。
例として、𝐴m と𝜎 2 の 2 種類のパラメータを利用し、2 次元格子点を形成した図を図1.5 に示す。今回求めるピンニングパラメータは 4 つであるため、4 次元の格子を設け
る。

図 0.5

2 つのパラメータを格子状に組み合わせた概略図

メッシュ法を用いる利点としては、パラメータの範囲を分割して計算するだけである
ためアルゴリズムの実装が簡単であり、全体を一様に計算するため、局所解に陥らないという利点がある。欠点としては、分割数が増大すると、計算時間が指数的に増大する欠点がある。

1.6

シンプソン法シンプソン法とは、ある関数𝑓(𝑥)をラグランジュ補間よって多項式に近似し、近似した多項式の積分値を求めることで、関数𝑓(𝑥)を数値的に積分する方法である。
𝑓(𝑥)の区間[𝑎, 𝑏]上の点を𝑛等分した場合、生成される多項式𝑃𝑛 (𝑥)は、- 10 -𝑛𝑃𝑛 (𝑥) = ∑ 𝑓(𝑥𝑗 )𝑙𝑗 (𝑥)(1.36)

𝑗=0(1.36)式の多項式𝑙𝑗 (𝑥)は、𝑙𝑗 (𝑥) =(𝑥 − 𝑥0 )(𝑥 − 𝑥1 ) … (𝑥 − 𝑥𝑗−1 )(𝑥 − 𝑥𝑗+1 ) … (𝑥 − 𝑥𝑛 )(𝑥𝑗 − 𝑥0 )(𝑥𝑗 − 𝑥1 ) … (𝑥𝑗 − 𝑥𝑗−1 )(𝑥𝑗 − 𝑥𝑗+1 ) … (𝑥𝑗 − 𝑥𝑛 )(1.37)

で表される。𝑛 = 3の場合の多項式𝑃3 (𝑥)を利用すると、積分値は、
𝑏𝑏∫ 𝑓(𝑥)d𝑥 ≅ ∫ 𝑃3 (𝑥)d𝑥
𝑎

=

𝑎

3ℎ
[𝑓(𝑎) + 3𝑓(𝑎 + ℎ) + 3𝑓(𝑎 + 2ℎ) + 𝑓(𝑏)]8(1.38)となる。ただし、ℎ = (𝑎 − 𝑏)/3とする。このときの概略を図 1.6 に示す。

図 0.6

積分区間を 3 等分した場合のシンプソン法

シンプソン法の誤差は、積分区間が広がるほど大きくなる。したがって、誤差を小さくするために、区間[𝑎, 𝑏]を3𝑡（𝑡は自然数）等分して、区間 (𝑥3𝑖 , 𝑥3𝑖+3 )(𝑥0 = 𝑎, 𝑥3𝑡 = 𝑏)にシンプソン法を適応し、その和を考えると、
𝑎𝑛−1∫ 𝑓(𝑥)d𝑥 ≅ ∑
𝑏

𝑖=0

3ℎ
{𝑓(𝑥3𝑖 ) + 3𝑓(𝑥3𝑖+1 ) + 3𝑓(𝑥3𝑖+2 ) + 𝑓(𝑥3𝑖+3 )}8- 11 -

𝑛−1

3(𝑎 − 𝑏)=[𝑓(𝑎) + 𝑓(𝑏) + ∑{3𝑓(𝑥3𝑖+1 ) + 3𝑓(𝑥3𝑖+2 ) + 2𝑓(𝑥3𝑖 )}]8𝑡(1.39)

𝑖=1のようになる。(1.39)式は係数が 3/8 となっているので、シンプソン 3/8 則と呼ばれる。

1.7

GPGPU (General-purpose computing on graphics processing units)

GPU(Graphics Processing Unit)はコンピュータに積まれているユニットで、主に画像処理や描写につかわれている。GPU の例として、本研究で使用する GeForce GTX 960 を
図 1.7 に示す。GPU には、画像処理に特化したプロセッサが多数詰まれており、一つ一つのプロセッサの性能は CPU に劣るが、大量にあるプロセッサを生かした並列処理が得意である。主に画像処理に使われている多数のプリプロセッサを計算やシミュレーシ
ョンなどに利用することを GPGPU と呼ぶ。

図 0.7

GeForce GTX 960

並列化できる箇所を GPGPU により並列化すれば、高速化は達成できる。しかし、
GPGPU による最適化は非常に難しい。GPGPU を扱う上で処理速度に影響しやすい箇所の例として、分岐処理、メモリアクセスの 2 つを挙げる。
GPU による分岐処理は、条件を考えずに、処理を実行し、条件が不適当な場合は入力や出力を行わないようにして正しい結果を保証している。しかし、処理が増えるので、処理時間が余分にかかってしまう。実際のプログラミングにおいては、分岐命令を減らすことや実行される頻度の低い部分はできるだけ短くするなどの工夫が必要である。
次に、データ転送について挙げる。GPU 内部で計算や処理を行う場合、CPU のメモ
リに GPU はアクセスできないため、データ転送を行う必要がある。そのため、大規模なデータが必要な場合にはデータ転送に相応の時間がかかってしまう。また、不連続な- 12 -

メモリアクセスを多くすると、それだけアクセス時間がかかってしまう。したがって、できるだけ連続したメモリ領域にアクセスするように変更すると処理時間の向上が見込める。また、GPU には条件を満たせばデータにまとめてアクセスできるコアレッシ
ングという機能が備わっている。コアレッシングを使うことにより、データ転送の効率化が見込める。このように、メモリアクセスについては考慮すべき点が多くある。

1.8

GPU の構造GPU を構成する部品としては、基盤、GPU チップ、DRAM(Dynamic Random Access
Memory)および画像出力端子、電源、冷却装置から構成されている。GPU において演算を担当している部分が GPU チップである。GPU の概略を図 1.8 に示す[4]。
GPU

SM (Streaming
Chip
SM

Multiprocessor)内部
SM

SM

L1 キャッシュ
SM

SM

SM

共有メモリ
L2 キャッシュ
テキスチャメモリ
コンスタントメモリ

グローバルメモリ

レジスタ
CUDA Core
図 0.8

GPU の構造

1.8.1 ストリーミングマルチプロセッサ（SM）
ストーミングマルチプロセッサとは、GPU チップ上に搭載されている演算部で、それぞれに、図 1.8 に示すとおりレジスタやキャッシュ、演算装置にあたる CUDA コア
がいくつも備え付けられている。CUDA コアとは、GPU における演算子の最小単位である[5]。

- 13 -

1.8.2 レジスタ
CUDA コアに最も高速にアクセスできるメモリである。主な用途として、カーネルで使用される変数のための領域である。スレッドごとに独立した領域を持ち、共有ができない。

1.8.3 共有メモリ
シェアードメモリとも呼ばれる、SM 内で共有されるメモリのことである。L1 キャッ
シュと容量を共有している。レジスタよりはメモリアクセス速度は劣るものの、アクセスが高速なメモリである。

1.8.4 テキスチャメモリ
テキスチャと呼ばれる最大 3 次元のデータを配置して、そのインデックス（添え字）を読み込むことでデータの読み込みが可能である。書き込みはできず、読み込み専用である。本研究のプログラムでは利用していないが、連続したデータを保持するために用いられることがある。

1.8.5 コンスタントメモリ
CPU 側からは書き込みができるが、GPU 側からは読み込みしかできないメモリである。コンスタントメモリの内容は全てコンスタントキャッシュに移動されるため、容量は非常に小さい。CUDA プログラミングにおいては、定数を保持するためによく用いられる。

1.8.6 グローバルメモリ
GPU の主メモリである。GPU メモリチップではなく、GPU の DRAM に備え付けられている。容量はレジスタや共有メモリと比較して非常に大きいが、低速である。

1.8.7 キャッシュ
SM 内部に備え付けられており、SM 内でデータが共有されるもの（L1 キャッシュ）と GPU 内部に備え付けられており、SM 間でデータが共有されるもの（L2 キャッシュ）がある。キャッシュは、メモリアクセス速度を改善するために使われる。

- 14 -

1.9

CUDA

CUDA とは、NVIDIA 社が提供している C 言語を拡張した総合開発環境である。主にコンパイラやライブラリなどから構成されている。他に GPGPU を利用する際に使われる総合開発環境としては、OpenCL や Direct Compute などがある。これらと比較すると、
CUDA は GPU の性能を引き出しやすくプログラミングできる。また、GPGPU を利用するための煩雑な準備を簡単な形で記述できるため、プログラマーはアプリケーションやアルゴリズムの実装に注力することができる。

1.10

GPU における並列化の階層

GPU による並列化の階層として、スレッド(Thread)、ブロック（Block）、グリッド（Grid）の 3 つの要素がある。スレッドはプログラムが動く最小の単位である。GPU においては、CUDA コアが該当する。
ブロックはスレッドを集めたものであり、1 ブロックには最大 512 のスレッドを集めることができる。ブロックは 3 次元まで拡張させて管理させることができ、1 ブロックに 512 スレッド、1 ブロックに 16×16 スレッド、1 ブロックに 8×8×8 スレッドのように管理することができる。
グリッドは、ブロックを集めたものである。ブロックと同様に 3 次元的に表すことができるが、𝑧方向は 1 に固定されている。グリッド数にも上限があり、𝑥方向と𝑦方向の積が 65536 となるように設定しなければならない。

1.11 実験の目的
超伝導体の応用のためには、超伝導体材料の電界-電流密度（𝐸-𝐽）特性を知ることが必要となる。𝐸-𝐽特性を求める手法の一つとして、磁束クリープ・フローモデルがある。
磁束クリープ・フローモデルの計算においては、最適なピンニングパラメータを求めることが必要であり、ピンニングパラメータを求める方法の一つとしてメッシュ法がある。
メッシュ法は、各パラメータを分割し、分割したパラメータを組み合わせて格子状にし、その格子上の全ての点について計算を行うので、局所解に陥らない利点があるが、パラメータの分割を細かくする程、計算時間が著しく増大する欠点がある。計算時間の削減には、並列処理が有効であることが知られている。並列処理の手法の 1 つとして
GPGPU がある。
本研究では、GPGPU でメッシュ法を計算することで計算の高速化および精度の向上を目的とする。

- 15 -

第2章 計算手法
2.1

評価値𝑑の計算方法

ピンニングパラメータの精度を表す評価値として、距離𝑑を用いる。距離𝑑は、(2.1)

log(𝐸exp ) − log(𝐸theo )𝑑=12∑(log(𝐸exp ) − log(𝐸theo ))
𝑀

(2.2)

とする。ここで、𝑀はサンプル数、𝐸exp は測定値、𝐸theo はピンパラメータを用いて計算した理論値である。図 2.1 に𝑑を模式的に表したものを示す。

図 2.1

距離𝒅の模式図

今回計算するピンニングパラメータは、𝐴m 、𝜎 2 、𝛾、𝑚の 4 つである。4 つのピンニ
ングパラメータとパラメータの分割数の範囲を表 2.1 に示す。
表 2.1

log𝐴m

パラメータの範囲

𝜎2

𝑚

分割数 N

−1

𝛾
−2

最大値

12

1.5 × 10

7.0 × 10

6.0

14

最小値

10

5.0 × 10−3

5.0 × 10−1

2.0

8

- 16 -

2.2

プログラムの流れ本研究で作成したアルゴリズムの概略を図 2.1 に示す。青色が CPU の動作、橙色がGPU の動作を示す。𝑁の部分のループは𝑚パラメータを設定して計算するためのものであり、𝑀の部分のループは、距離𝑑を求めるためのループである。

GPU に利用するデータ領域を確保

ピンニングパラメータ設定

GPU のデータをコピー

𝑖 = 1; 𝑖 < 𝑁; 𝑖 + +終了

𝑗 = 1; 𝑖 < 𝑀; 𝑗 + +式 2.1 式を計算

𝑖が𝑀になるまで繰り返す

𝑖が𝑁になるまで繰り返す

図 2.2

プログラムの流れ- 17 -

2.3

実行環境

実行環境を表 2.2 に示す。
表.2.2

2.4

実行環境

CPU

Intel(R) Core(TM) i5-6500 CPU

GPU

NVIDIA GeForceGTX 960

解析に用いるデータ

今回解析に用いた線材の諸元を表 2.3 に示す。この線材は、国際超電導産業技術研究センター(ISTEC/SRL)によって作製された線材で、GdBa2 Cu3 O7−𝛿 (GdBCO)コー
ト線材にBXOピンを導入したものである。温度は、20 K、30 K のそれぞれの温度の条件下で磁束密度が 1 T、2 T、3 T の 3 種類のデータを用いて、CPU、および GPU によって計算を行う。
表.2.3

使用する線材の諸元

臨界温度 𝑇c ⁡ [K]90.3

0 K における上部臨界磁界 𝐵𝑐2 ⁡ (0⁡ K)[T]100

常伝導状態における抵抗率 𝜌𝑛 [Ω ∙ m]

1.5 × 10−6 ⁡

ピ

6.28

ン
超

1.0 × 10−6 ⁡

形
伝
状
導
に
層
依
の
る
厚
𝑡⁡ 定
[m]
𝜁

- 18 -

第3章 結果・考察
3.1

CPU と GPU の計算時間の比較

図 3.1 では、縦軸には、秒を単位とした計算時間を表しており、横軸にはパラメータの分割数を表している。丸のシンボルが CPU の結果であり、四角のシンボルが GPU の結果である。CPU の計算時間の近似式を赤線で、GPU の計算時間の近似式を青線で示している。近似式はパラメータが 8 から 14 までの実行時間で近似したものである。
分割数が 8 の場合には CPU と GPU の計算時間は 200 秒程度と近い値をとっている。
しかし、分割数が 14 の場合には、CPU の計算時間がおよそ 1700 秒であるが、GPU の計算時間はおよそ 260 秒となっている。したがって、CPU の計算時間の増大は、GPUの計算時間の増大より大きい。分割数が 8 から 14 における CPU と GPU の計算時間と分割数の関係を近似式で表すと、式 4.1、式 4.2 のようになる。

1000
C
PU

実行時間 [s]

1500

500
GPU

0

5

10

パラメータの分割数 N
図 3.1

パラメータの分割数と CPU、GPU の計算時間- 19 -

CPU: 4.00 × 10−2 𝑁 4.04

(3.1)

GPU: 1.39 × 10𝑁1.22(3.2)

近似式から、CPU も GPU も計算時間は指数関数的に増大していることが分かる。しかし、𝑁の次数を比較すると、GPU で計算した場合の近似式の次数が小さい値となっているため、GPU で効率よく計算できることが分かる。分割数をさらに増やした場合、結果から GPU と CPU の計算時間の差はさらに広がると予想されるので、GPU での計算
は分割数が増大するほどに有効であると考えられる。CPU の近似式の係数は GPU の近似式の係数より小さいから、分割数が小さい場合は CPU の計算が GPU より早いと考えられる。
パラメータの分割数と評価値𝑑の関係を図 3.2 に示す。
縦軸は評価値であり、横軸は 8 から 14 までの分割数である。分割数が 8 から 14 まで
増加する間に、評価値は、8.36 × 10−2 から3.81 × 10−2 まで減少している。全体の傾向としては、分割数が増大するほど、評価値は小さくなる。したがって、分割数が増大すると、ピンニングパラメータを精度よく求めることができる。
分割数が 9 の時だけ、全体の傾向から外れている。これは、分割数が小さい場合でも、

評価値d

0.1

0.05

0
8

10

12

メッシュ法における分割数 N
図 3.2

パラメータの分割数と評価値の関係- 20 -

14適切なピンニングパラメータに近いパラメータの組み合わせがあれば評価値は小さくなることを示している。
分割数が 8、11、14 のそれぞれについて最適なピンニングパラメータを用いて求めた
𝐸-𝐽特性と実験の特性を比較したものを図 3.3 に示す。また、その時の𝑑の値を表 3.1 に
示す。

1T
2T
3T

電界 E [V/m]

10

10

−9

30 K

−10

20 K

分割数8
分割数11
分割数14

10

−11

10

2

電流密度 J [A/m ]
図 3.3

11

分割数を変化させた場合の磁束クリープ・フローモデルによる理論値の𝑬-𝑱特性と実験値の𝑬-𝑱特性との比較
表 3.1

分割数
評価値𝑑

分割数が 8、11、14 の時の𝒅の値

8118.36 × 10−2

5.84 × 10

14−2

3.81 × 10−2

図 3.3 では、縦軸に電界をとり、横軸に電流密度をとったグラフとなっている。菱形のシンボル、丸のシンボル、正方形のシンボルがそれぞれ 1 T、2 T、3 T における実験値となっており、20 K、および 30 K の条件における実験結果は線分で纏めることで示している。
赤線は分割数が 8 の場合の最適なピンニングパラメータで計算した𝐸-𝐽特性、青線は- 21 -

分割数が 11 の場合の最適なピンニングパラメータで計算した𝐸-𝐽特性、緑線は分割数が14 の場合の最適なピンニングパラメータで計算した𝐸-𝐽特性である。
分割数が 8 の場合、1 T、20 K の条件における実験値と計算した𝐸-𝐽特性を比較すると差が生じている。しかし、分割数が 11、分割数が 14 となるにつれ、実験値と計算した
𝐸 -𝐽特性の差が小さくなっている。別の温度、磁束密度の条件下においても、分割数を多くして推定したピンニングパラメータで計算した𝐸 -𝐽特性と実験値との差が小さくなる傾向がある。したがって、分割数が多いほど実験値に合ったピンニングパラメータが求められていることが分かる。

3.2

GPU で一度に計算する条件について

GPU の計算において、{𝐴𝑚 , 𝜎 2 , 𝛾}のパラメータを並列化しており、𝑚は、パラメータを代入する際に引数として決定している。ピンニングパラメータから電界を求めるには、他に𝐵, 𝑇, 𝐽のパラメータが必要である。3.1 節において、𝐵, 𝑇, 𝐽のパラメータは GPU
で計算を行う際には、引数として渡している。
電界を求める計算においては、𝐵, 𝑇のパラメータとピンニングパラメータだけで計算できる箇所が存在する。そのため、𝐽のパラメータを GPU に保持して一度に計算を行うことによって、計算が行えるのではないかと考えた。これが prog2 である。
prog2 では、𝐽のデータを GPU で保持しているが、CPU で処理を戻す度に、𝐽のデータ
を CPU から GPU にコピーしている。したがって、始めに全ての𝐽のデータを GPU にコ
ピーし、𝐵, 𝑇の条件も GPU に保持して計算することによって、prog2 から計算の効率化を図ったプログラムが prog3 である。GPU 内部に保持するパラメータを可変のパラメ
ータ、GPU 内部で引数として与えて 1 度の GPU の計算内では変わらないパラメータを固定するパラメータと呼ぶこととする。これらを纏めたものを表の 3.2 に示す。
表 3.2

それぞれのプログラムにおける計算条件

固定するパラメータ

可変のパラメータ

Prog1

𝐵, 𝑇, 𝐽

なし

Prog2

𝐵, 𝑇

𝐽

Prog3

なし

𝐵, 𝑇, 𝐽- 22 -

1000
prog1
prog2

計算時間 [s]

prog3

500

0
8

10

12

14

パラメータの分割数 N
図 3.4

GPU で一度に計算する条件を変更した結果

表 3.2 の条件の下での、計算結果を図 3.4 に示す。縦軸は単位を秒とした実行時間を表しており、横軸は 8 から 14 までのパラメータの分割数を表している。青のユニット
が prog1 の実行時間を表しており、青の実線は prog1 の計算時間の近似式を表している。緑のユニットは prog2 の実行結果を表しており、緑の実線は prog2 の計算時間の近似式を表している。赤のユニットは prog3 の実行結果を表しており、赤の実線は
prog3 の計算結果の近似式を表している。
求めた近似式を、式 4.3、式 4.4、式 4.5 で示す。

prog1:1.39 × 10𝑁1.22

4.3

prog2:1.29 × 10𝑁1.21
prog3:2.23 × 10𝑁1.11

4.4
4.5

prog2 の計算時間は、全ての分割数において実行時間は一番短く、prog3 の計算時間が一番長い結果となっている。prog1 では、𝐽 の値によらない箇所も毎回計算していたが、prog2 では、𝐽 の値によって変化しない箇所は 1 度しか計算していないので、
prog2 の計算回数は、prog1 の計算回数より少ない。計算時間が少なくなったことによって、prog2 の計算時間は、prog1 の計算時間より短くなったと考えられる。しかし、
prog2 と同じように計算回数を減らした prog3 では逆に計算時間が長くなった。prog3
では GPU で動作する時間が長いため、その分 GPU に負荷がかかり、計算時間に差が- 23 -

出てしまったと考えられる。他に考えられる原因としては、prog3 のプログラムはルー
プ回数が prog1、prog2 よりも多い。ループ回数以外のプログラムには差異はないため、ループ処理に時間がかかって、計算時間が増加したと考えられる。

- 24 -

第4章 まとめ磁束クリープ・フローモデルは、𝐸-𝐽特性を理論的に求める際に利用されている。磁
束クリープ・フローモデルを使うにあたって、最適なピンニングパラメータの推定が必要となる。
ピンニングパラメータの推定にはメッシュ法が用いられている。ただし、メッシュ法
では、分割数を増やしていくと精度は向上するが、計算時間は著しく増大する欠点がある。その欠点を解消するため、GPGPU による計算の並列化を行う。
GPGPU とは、普段は画像処理に使われるプリプロセッサを画像処理以外の目的に利用することである。多数のプリプロセッサが計算やシミュレーションの並列化を可能にしている。
GPGPU でメッシュ法を計算することで計算の高速化および精度の向上を目的として、研究を行った。
CPU と GPU も分割数の増大によって、計算時間は増大している。ただし、パラメー
タの分割数と CPU、GPU の計算時間の関係を近似式で表すと、分割数の増大とともにGPU の計算時間の増加は CPU の計算時間の増加と比較すると小さいことがわかる。したがって、GPU による計算の高速化が確認できた。また、分割数を増やしていくほど、GPU と CPU の計算時間の差が広がる。分割数を増やして計算を行う場合は、GPU を用
いて計算を行ったほうが、計算時間を短縮できる可能性がある。
分割数を増加させると、評価値𝑑は小さくなる傾向にある。よって、パラメータの精度を向上させるためには、分割数を増加させていくことが必要であることが分かる。
磁束クリープ・フローモデルにおいては、GPU で一度に計算する条件を変更することによって、計算の一部を省略でき、高速化できる。また、GPU で一度に処理する時間が長い場合やループ処理が多くなる場合は、計算時間が増大する可能性がある。

- 25 -

謝辞
本研究を進める際に、熱心なご指導や適切な助言をしていただいた九州工業大学情報工学部電子情報工学科小田部荘司教授、木内勝准教授に心より感謝致しております。
小田部教授には、常日頃から研究者としてのあり方や、研究に関する心構えやルールなどを細かくご指導いただきました。深く感謝申し上げます。
木内准教授には、超伝導に関する知識や研究者としてのあり方をご教授していただきました。深く御礼を申し上げます。
また、九州工業大学大学院先端情報工学専攻電子情報工学分野の田邊裕也氏には、論文作成から平時の資料作成まで、数多くのご助力を頂きました。心から感謝を申し上げます。
最後に、公私共々お世話になりました、小田部・木内研究室の皆様に深く御礼申し上げます。

- 26 -

参考文献[1]松下照男(1994)「磁束ピンニングと電磁現象」産業図書[2] 増田嘉道(2013)「最急降下法を用いた超伝導体の電界-電流密度特性評価の検討」九
州工業大学情報工学部電子情報工学科卒業論文[3]原田將敬(2016)「超伝導体の電界‐電流密度特性の評価における計算最適化」九州工業大学情報工学部電子情報工学科卒業論文[4]青木尊之、額田彰(2009)「はじめての CUDA プログラミング」工学社[5]秋吉祐弥(2012)「遺伝的アルゴリズムを用いた超伝導体の電界-電流密度特性評価の
GPGPU に よる高速化に関する研究」九州工業大学情報工学部電子情報工学科卒業論
文- 27 -