TDGL方程式に基づく数値シミュレーション
オブジェクト指向プログラミングによる込み入った数式の実装（前回とほぼ同じ内容でごめんなさい．．．）
有明高専 電子情報工学科 松野哲也
「量子化磁束動力学シミュレーション研究会 」（2016年１月１９日（火）14:00-16:00）
九州工業大学サテライトキャンパス（天神イムズ11F）

1

FDTD法によるTDGLシミュレーション● FDTD： 時間領域有限差分法（ 偏微分方程式 → 時間発展差分方程式 ）● FEM ： 有限要素法（ 偏微分方程式 → 変分問題 → 線形方程式 ）2

TDGL方程式の導出

3

TDGL方程式
時間依存ギンツブルグランダウ方程式( Time dependent GinzburgLandau equation )4

ゲージ対称性5

なぜ φゼロゲージ，なぜリンク変数法か？

FDTD法と相性が良い．

flux flow シミュレーションを行うにはリンク変数法が必須．

陽的数値積分を組める．
並列化が容易高速化

E 一定は A の時間微分が一定であることを意味するから．
つまり，t → ∞ で A は有界ではない．

6リンク変数法
Ny

格子点の上でオーダーパラメータを定義.

格子点と格子点の間でベクトルポテンシャルを定義.
.

2112. . .

Nx
7yUxi, j+1i, j+1Lz i, j

Uy i, ji, j

Ux i, j

リンク変数：

Uy i+1, j

hxi+1, jh
Ny + 1
Ny.
.
.

2
100

12. . .

Nx

Nx + 1

ψUx

Uy8

共変微分演算子（２次元）

9

Uxi, j+1
リンク変数の積は経路積分に対応する．

Uy i, ji, j

Uy i+1, j

Ux i, j

10z

3次元では

Ly i, j, kx

Lx i, j, ki, j, kyLz i, j, k

11

２次元

Lz i-1, j

Lz i, jLz i, ji, j

Lz i, j-1i, j

12Semi-discretized equation ( 2次元 )13

Semi-discretized equation ( 3次元 )14

コーディング（実装）
表現に類似性をもたせる
a = a + b・c・d + e・fa.equal( a.plus( b.mul(c).mul(d) ).plus( e.mul(f) );Complex plus( Complex z ){return new Complex( x + z.x, y + z.y );}15

「複素数クラス」を定義

Complex mul(Complex z){
double x1 = z.x;
double y1 = z.y;
double x2 = x;
double y2 = y;return new Complex(x1*x2 - y1*y2, x1*y2 + y1*x2);}16

Complex Lz(Complex[][][] Ux, Complex[][][] Uy, Complex[][][] Uz, int i, int j, int k){
return Ux[i][j][k].mul(Uy[i+1][j][k]).mul(Ux[i][j+1][k].conj()).mul(Uy[i][j][k].conj());}

17

//////////
void do_Laplacian(Complex[][][] psi, Complex[][][] Ux, Complex[][][] Uy, Complex[][][] Uz, double dt, double D){for(int k = 1; k <= Nz; k++){for(int j = 1; j <= Ny; j++){for(int i = 1; i <= Nx; i++){psi_s[i][j][k].equal(psi[i][j][k].plus(((psi[i+1][j][k].mul(Ux[i][j][k])).plus(psi[i-1][j][k].mul(Ux[i-1][j][k].conj())).plus(psi[i][j+1][k].mul(Uy[i][j][k])).plus(psi[i][j-1][k].mul(Uy[i][j-1][k].conj())).plus(psi[i][j][k+1].mul(Uz[i][j][k])).plus(psi[i][j][k-1].mul(Uz[i][j][k-1].conj())).minus(psi[i][j][k].mul(6.0) ) ).mul(dt * D / (h * h) )));
ラプラシアン演算の実装}}}
for(int k = 1; k <= Nz; k++){for(int j = 1; j <= Ny; j++){for(int i = 1; i <= Nx; i++){
psi[i][j][k].equal(psi_s[i][j][k]);}}}}//////////
18

境界条件
境界での磁束密度指定
(b)

(a)

hj = jB + 1

j = jB + 1hj = jB
i=0

i=1

j = jB

i=2(c)

i = Nx - 1 i = Nx i = Nx + 1(d)j = Ny + 1i = iB

j = Nyj=2
j = Ny - 1j=1i = iB

j=019横磁界

B

縦磁界

B, JB, JJBJ20境界条件の実装/// constraints by the applied field and current/// y-z plane ( Bz-Uy-component ) set Ja = 0 when longitudinalfor(int k = 1; k <= Nz; k++){for(int j = 1; j <= Ny-1; j++){ex.setPolar(1.0, h*h*(Ba + Ja*SSx/2.0));
Uy[1][j][k].equal(Uy[2][j][k].mul(Ux[1][j+1][k].conj()).mul(Ux[1][j][k]).mul(ex));ex.setPolar(1.0, -h*h*(Ba - Ja*SSx/2.0));
Uy[Nx][j][k].equal(Uy[Nx-1][j][k].mul(Ux[Nx-1][j+1][k]).mul(Ux[Nx-1][j][k].conj()).mul(ex));}}/// z-x plane ( Bz-Ux-component ) omit when transverse
for(int i = 1; i <= Nx-1; i++){for(int k = 1; k <= Nz; k++){ex.setPolar(1.0, - h*h*Ba);Ux[i][1][k].equal( (Uy[i+1][1][k].conj()).mul(Ux[i][2][k]).mul(Uy[i][1][k]).mul(ex) );ex.setPolar(1.0, h*h*Ba);Ux[i][Ny][k].equal( Ux[1][Ny-1][k].mul(Uy[i+1][Ny-1][k]).mul(Uy[i][Ny-1][k].conj()).mul(ex) );}}/// y-z plane ( J=(0,0,Jb)-Uz-component ) longitudinalfor(int k = 1; k <= Nz-1; k++){for(int j = 1; j <= Ny; j++){ex.setPolar(1.0, h*h*Jb*SSx/2.0);
Uz[1][j][k].equal( Ux[1][j][k].mul(Uz[2][j][k]).mul(Ux[1][j][k+1].conj()).mul(ex) );Uz[Nx][j][k].equal( (Ux[Nx-1][j][k].conj()).mul(Uz[Nx-1][j][k]).mul(Ux[Nx-1][j][k+1]).mul(ex) );}}/// z-x plane ( J=(0,0,Jb)-Uz-component ) longitudinalfor(int i = 1; i <= Nx; i++){for(int k = 1; k <= Nz-1; k++){ex.setPolar(1.0, h*h*Jb*SSy/2.0);
Uz[i][1][k].equal( Uy[i][1][k].mul(Uz[i][2][k]).mul(Uy[i][1][k+1].conj()).mul(ex) );Uz[i][Ny][k].equal( (Uy[i][Ny-1][k].conj()).mul(Uz[i][Ny-1][k]).mul(Uy[i][Ny-1][k]).mul(ex) );}}21

シミュレーション例
••
•
•
•

FEM（φ given）（p00）2次元横磁界（φ given）（p01）2次元横磁界（φ given）（p02）2次元横磁界（リンク変数，φゼロゲージ）（p03）＊電場表示
3次元縦磁界（リンク変数，φゼロゲージ）（p04）＊電場の表示はリンク変数の偏角の時間微分を利用

22

エネルギーEK + EP + EB

EK + EP
EBEK
EPtime23

今後の課題・展開• アルゴリズムの検討・改良• 並列化（高速化）e.g. by GPGPU ( Java + OpenCL )• ピンニング• Two-component 系

24