FDTD法によるTDGLシミュレーション
リンク変数法によるシミュレーション，および オブジェクト指向プログラミングによる込み入った数式の実装
有明高専 電子情報工学科
松野哲也
応用物理学会２０１５秋（９月１３日（日）12:00-­‐13:30）　インフォーマルミーティング　「量子化磁束動力学シミュレーション研究会 」
1
FDTD法によるTDGLシミュレーション
● FDTD： 時間領域有限差分法 （ 偏微分方程式 → 時間発展差分方程式 ） ● FEM ： 有限要素法 （ 偏微分方程式 → 変分問題 → 線形方程式 ）
2
ix A
Time dependent GL equation
TDGL方程式の導出
ensity is defined as
ϵGL = ϵK + ϵP + ϵB ,
ϵK
ϵP
ϵB
1
=
|Dψ|2 , D ≡ ∇ − igA,
2
η
=
(|ψ|2 − 1)2 ,
2
1
=
|∇ × A|2 .
2
iation:
δ
≡
∂
−∇·
∂
,
δ
=
∂
−∇·
∂
3
, i = x, y, z
∂t
+ ∇φ = −
δA
,
e TDGL equations will be obtained as
!
" us note that the equations (25) and (
alar potential.
Let
∂
δϵGL
τ
+ igφ ψ = −
,
(25
transformation:
∂t
δψ
TDGL方程式 ∂A
δϵGL ∂χ 時間依存ギンツブルグランダウ + ∇φ φ=−→
− φ −,
(26
方程式 A −→ A∂t
+ ∇χ,
,
ψ
−→
ψ
exp(igχ).
δA
( Time dependent Ginzburg-­‐
∂t
)
tential. Let us note that the equations (25)Landau andequa@on (26) are
invarian
Lrmation:
equations is obtained as follows:
!
"
∂
∂χ
→A+
φ −→
− =,
τ ∇χ, +
igφ φ ψ
∂t
∂t
∂A as follows:
ions is obtained
+ ∇φ =
"
∂t1
+ igφ ψ =
2
1 2
2
ψD
−→ψψ+
exp(igχ).
(27
η(1 − |ψ| )ψ,
2
#
$
gIm ψDψ − ∇ × ∇ × A.
(∇ − igA)2 ψ + η(1 − |ψ|2 )ψ,
#
$
B+ ∇φBoundary
condition
= gIm ψ(∇ − igA)ψ − ∇ × ∇ × A.
A
(28
4
(29
∂t
+ ∇φ = −
δA
,
he scalar potential. Let us note that the equations (25) and (26)
auge transformation:
ゲージ対称性
equation
expressed
byvariables
link varia
n expressed
by link
∂χ
A −→ A + ∇χ, φ −→ φ −
, ψ −→ ψ exp(igχ).
∂t
TDGL
equations
is obtained
asgauge
follows:is as follows:
equations
of
the
φ-zero
the φ-zero
! gauge
" is as follows:
∂
1 2
2
τ ∂ψ + igφ 1
ψ =2 D ψ + η(1 − |ψ|
)ψ,
2
1τ ∂t2 = D ψ2+2η(1 − |ψ| )ψ,
=
D∂tψ
#
$
2 − |ψ| )ψ,
∂A+ η(1
− ∇ × ∇ × A.
2 ∂A ∂t + ∇φ =! gIm ψDψ
"
! = gIm
"
ψDψ − ∇ × ∇ × A.
∂t ψDψ − ∇ × ∇ × A.
= gIm
dix B Boundary condition
y discretized version of the above eqs.(14) and (15)
a rectangle
the two-dimensional
simulation
area. There
are l
dderversion
ofasthe
above eqs.(14)
and (15)
can be
r parameter
5
なぜ φゼロゲージ，なぜリンク変数法か？
on expressed
by link variables
of the φ-zero gauge is as follows:
ψ
1FDTD法と相性が良い．
2
2flux flow シミュレーションを =
D ψ + η(1 − |ψ| 行うにはリンク変数法が必須．
)ψ,
t
2
!
"
A
陽的数値積分を組める．
= gIm
ψDψ − ∇ × ∇ × A.
t
E 一定は A の時間微分が一定 であることを意味するから． つまり，t → ∞ で A は有界ではない．
並列化が容易
ed version 高速化 of the above eqs.(14) and (15) can be
i+1,j
i−1,j
i−1,j
i,j
i,j+1
i,j−1
i,j−16
リンク変数法
Ny
格子点の上で オーダーパラメータを定義
.
格子点と格子点の間で ベクトルポテンシャルを定義
.
.
2
1
1
2
. . .
Nx
7
Appendix C
Equations with the link v
y
j+1 link variables as
Let us defineUxi,the
リンク変数：
i, j+1
Lz i, j
Uy i, j
i, j
Ux i, j
Uy i+1, j
Uxi,j,k = exp(−ihgAi,j,k
x ),
Uyi,j,k = exp(−ihgAi,j,k
y ),
h
x
i+1, j
Uzi,j,k = exp(−ihgAi,j,k
),
z
h
where h is the
lattice spacing. By using these link variables, w
discrete version of the equations (28) and (29).
The dynamics of the order parameters are described as follows
Ny + 1
Ny
.
∂ψ i,j,k
1
i−1,j,k i−1,j,k
i,j,k i+1,j,k
i,j,k i,j+1
τ
=
(U
ψ
+
U
ψ
+
U
ψ
x
x
y
2
∂t
2h
i,j,k−1 i,j,k−1
+Uzi,j,k ψ i,j,k+1 + U z
ψ
− 6ψ i,j,k ) + η(1
.
.
2
1
0
0
1
2
. . .
Nx
Nx + 1
ψ
Ux
8
Uy
8
n
1
Dx ψi,j=1 i+1,j
−
igA
= Vi,jx ni,j(Vi,j
i,j
xi+1,jψi,j
x ψ),
i+1,j
1
i,j
i+1,j
i+1,j
i+1,j
i,jψ i,j )
ψ
−→
V
(V
ψ
−
V
∂x
∂x
−→D
V
(V
ψ
−
V
ψ
)
x
Dx ψx −→x Vxx h (Vxx
ψ
− Vxx ψ )
x
h
h
1
i,j
1
i,j
i+1,j
i+1,j
i,j
1
共変微分演算子（２次元）
i,j
i+1,j
i+1,j
i,j
his,=we can(Vobtain
the
discretized
version
of
the
i+1,j
i+1,j
i,j
=
(V
V
ψ
−
ψ
)
ψxx Vxx− ψψ ) − ψ )
x=Vx h (V
h
h
1 i,ji+1,j
11 i,j i,j
1 i,j i+1,j
i+1,j i+1,j
i,j i,j
i,j
i,j
i+1,j
i,j
D
ψ
−→
V
(V
ψ
−
V
)
=
(U
ψ
−
ψ
),
= x (Ux=ψ
x−
x ψ
xxψψx),
(U
−
ψ
),
hh h
h
1
1
i,j i+1,j i+1,j
i,j
i,j
22 i,j 1
i−1,j
i−
1
i,j
i+1,j
i+1,j
i−1,j
i−1,j
i,
i+1,j
i+1,j
i−1,j
i−
=
(V
V
ψ
−
ψ
)
D
ψ
−→
V
(V
+
V
ψ
−→DxxVψx −→
(Vx Vxx ψ2x (Vxx+ Vxψ ψ + Vxx− 2Vψx
2
h hh2
h
111 i,j i−1,j
i−1,j
1 i,j i+1,j
i+1,j
i,j
i−1,j
i−1,j
i−1,j
i,j
i,j
i+1,j
i−1,j
=
(U
ψ
−
ψ
),
=
+
U
ψ
−
2
=
(U
ψ
+
U
ψ
−
2ψ
)
x
x
=
(U
ψ
+
U
ψ
−
x
x
22
2
x
x
h
h
h
h
i,j 1
2
i+1,j i+1,j
i−1,j i
D
ψ
−→
V
(V
ψ
+
V
ψ
plains
the
above
discretizations.
x
x
x
x
e
above
discretizations.
xplains the above discretizations.
h2
1 i,j i+1,j
i−1,j i−1,j
9
Uxi, j+1
リンク変数の積は経路積分に対応する．
Uy i, j
Uy i+1, j
i, j
define a variable Li,j
using
link
variables
as
z
Ux i, j
i,j
z
Li,j
z
≡
=
≃
i,j i+1,j i,j+1 i,j
Ux Uy U x U y
i,j
exp[−ihg(Ai+1,j
−
A
y
y
exp[−ih2 g(∇ × A)i,j
z ],
i,j
− Ai,j+1
+
A
x
x )]
is the z-component of the rotation of the vector potenti
10
) or (xi + h/2, yj + h/2). That is, the variable Li,j
correspo
z
α=1,2
1 i,j,k i−1,j,k i,j,k i,j,k−1
i, j, k
− 2 (Lz Lz
Lx Lz z
−Lx1),
h
i,j,k3次元では
!
∂Uz
i,j,k i,j,k i,j,k+1
2 i,j,k 2
= −ig
(1 + κ |ψα | )Im[ψ α Uz ψα
]
∂t
α=1,2
1 i,j,k i,j−1,k i,j,k i−1,j,k
− 2 (Lx LxLy i, j, k Ly Ly i, j, k − 1),
h
y
i,j,k
i,j,k x
p” variables Li,j,k
,
L
,
L
are defined as
x
y
z
Lz i, j, k
Li,j,k
x
=
Li,j,k
=
y
Li,j,k
=
z
i,j,k i,j+1,k i,j,k+1 i,j,k
Uy Uz
Uy
Uz ,
i,j,k i,j,k+1 i+1,j,k i,j,k
Uz Ux
Uz
Ux ,
i,j,k i+1,j,k i,j+1,k i,j,k
Ux Uy
Ux
Uy .
physical meaning of the loop variable is expponential of the
pa
11
i,j
i,j+1
i,j
= exp[−ihg(Ai+1,j
−
A
−
A
+
A
y
y
x
x )]
2
i,j
≃
exp[−ih
g(∇
×
A)
],
z
２次元
Lz i-1, j
Lz i, j
s the z-component of the
Lz i, jrotation of the vector po
r (xi + h/2, yj + h/2). That is, the variable Li,j
z corr
j-1
i, j
(xi + h/2, yj +i, jh/2). See
Lz i, Fig.3(a).
an obtain the expressions as
i,j i,j−1
Lz Lz
i,j i−1,j
Lz Lz
3
≃ exp[−ih g(∇ × ∇ ×
≃ exp[−ih3 g(∇ × ∇ ×
i,j
A)x ],
A)i,j
y ],
(12) are defined at the sites (xi +h/2, yj ) and (xi , yj
e same positions for Uxi,j and Uyi,j , respectively.
12
L equation expressed by link variables
GL equations of the φ-zero gauge is as follows:
2
TDGL equation expressed by link variables
∂ψ
1
∂t
2
∂ψ ! 1 2"
∂A
τ
=
D ψ + η(1 − |ψ|2 )ψ,
= gIm
ψDψ
− ∇ × ∇ × A.
∂t
2
∂t
!
"
∂A
2 φ-zero gauge is 2as follows:
The set of TDGL
τ equations
= ofDthe
ψ + η(1 − |ψ| )ψ,
∂t
(
(14
(
= gIm ψDψ − ∇ × ∇ × A.
(15
tially
discretized version of the above eqs.(14) and (15) can be obtained
Semi-­‐discre@zed equa@on ( 2次元 )
Hence, the spatially discretized version of the above eqs.(14) and (15) can be obtained as
1 i,j i,j i+1,j
i−1,j i−1,j
i,j−1 i,j−1
i,j i,j+1
i,j
∂ψ
1
i−1,j
i,j−1
= τ 2 (Ux =ψ
+
U
ψ
+
U
ψ
+
U
ψ
−
4ψ
i,j i+1,j
i−1,jy
i,j i,j+1 y
i,j−1
i,j )
x
(U
ψ
+
U
ψ
+
U
ψ
+
U
ψ
−
4ψ
)
x
x
y
y
2h ∂t
2
2h
i,j 2
+η(1 − |ψ +η(1
| )ψ−i,j|ψ, i,j |2 )ψi,j ,
(16(
∂Uxi,j
1 i,j−1 i,j
i,j i,j i+1,j 1i,j
i,j
i,j−1
i,j
i+1,j
i,j
i,j Lz − 1)U
i,jxi,j ,
= −igIm[ψ
U]U
− 2 (LL
(17(
x ψ − ]Ux (L
z
= −igIm[ψ
U
ψ
−
1)U
,
x
x
z
z
x
∂t
h
2
h
i,j
∂Uy
1 i,j i−1,j
i,j i,j i,j+1
i,j
i,j
=i,j −igIm[ψ
U
ψ
]U
−
(L
L
−
1)U
,
(18
1
i,j
y
y
z
z
y
2 i−1,j
i,j i,j+1
i,j
i,j
∂t
h
= −igIm[ψ Uy ψ
]Uy − 2 (Lz Lz
− 1)Uy ,
(
h
y using link variables and
ariables and
∂Uxi,j
∂Ai,j
x
i,j−1
i,j
i,j
Lz
Li,j
z −1
13
x
Uyi,j,k
Uzi,j,k
=
=
x
exp(−ihgAi,j,k
y ),
exp(−ihgAi,j,k
),
z
(43)
(44)
equa@on ( 3次元 hereSemi-­‐discre@zed h is the lattice
spacing.
By)using these link variables, we can get the spatially
iscrete version of the equations (28) and (29).
The dynamics of the order parameters are described as follows:
∂ψ i,j,k
1
i−1,j,k i−1,j,k
i,j−1,k i,j−1,k
i,j,k i+1,j,k
i,j,k i,j+1,k
τ
=
(U
ψ
+
U
ψ
+
U
ψ
+
U
ψ
x
y
y
∂t
2h2 x
i,j,k−1 i,j,k−1
i,j,k i,j,k+1
i,j,k
+U
ψ
+
U
− 6ψ
) + η(1 − |ψ i,j,k |2 )ψ i,j,k .
(45)
The dynamics of thez link variablesz are ψ
described
as follows:
∂Uxi,j,k
1 i,j,k i,j,k−1 i,j,k i,j−1,k
i,j,k i,j,k i+1,j,k
= −igIm[ψ Ux ψ
Lz Lz
− 1),
8] − 2 (Ly Ly
∂t
h
∂Uyi,j,k
1 i,j,k
i,j,k
i,j,k i,j,k−1
= −igIm[ψ Uyi,j,k ψ i,j+1,k ] − 2 (Lz Li−1,j,k
L
− 1),
z
x Lz
∂t
h
∂Uzi,j,k
1 i,j,k
i,j,k
i,j,k i−1,j,k
= −igIm[ψ Uzi,j,k ψ i,j,k+1 ] − 2 (Lx Li,j−1,k
L
− 1),
x
y Ly
∂t
h
(4
(4
(4
i,j,k
i,j,k
here the ”loop” variables Li,j,k
are defined as
x , Ly , Lz
i,j,k+1
Li,j,k
= Uyi,j,k Uzi,j+1,k U y
x
Li,j,k
y
=
i,j,k
Uz
,
i+1,j,k i,j,k
Uzi,j,k Uxi,j,k+1 U z
Ux ,
(4
14
(5
境界条件
境界での磁束密度指定
(b)
(a)
h
j = jB + 1
j = jB + 1
h
j = jB
i=0
i=1
j = jB
i=2
(c)
i = Nx - 1 i = Nx i = Nx + 1
(d)
j = Ny + 1
i = iB
j = Ny
j=2
j = Ny - 1
j=1
i = iB
j=0
15
数値シミュレーション
•  動画：２次元 横磁界 •  動画：３次元 横磁界，縦磁界，リング •  リアルタイムデモ ２次元(40x80)，３次元(30x30x30) ※ order parameter 値に応じた色で可視化 16
コーディング（実装）
17
表現に類似性をもたせる
a = a + b・c・d + e・f a.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 ); } 18
「複素数クラス」を定義
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); } 19
i,j,k
i,j,k
p” variables Li,j,k
,
L
,
L
are defined as
x
y
z
Li,j,k
x
=
Li,j,k
=
y
Li,j,k
=
z
i,j,k i,j+1,k i,j,k+1 i,j,k
Uy Uz
Uy
Uz ,
i,j,k i,j,k+1 i+1,j,k i,j,k
Uz Ux
Uz
Ux ,
i,j,k i+1,j,k i,j+1,k i,j,k
Ux Uy
Ux
Uy .
physical meaning of the loop variable is expponential of the pa
Complex Lz(Complex[][][] Ux, Complex[][][] Uy, Complex[][][] Uz, int i, int j, int k){ ntial.
return Ux[i][j][k].mul(Uy[i+1][j][k]).mul(Ux[i][j+1][k].conj()).mul(Uy[i][j][k].conj()); } 3
20
The dynamics of the order parameters are described as follows:
∂ψ i,j,k
1
i−1,j,k i−1,j,k
i,j−1,k i,j−1,k
i,j,k i+1,j,k
i,j,k i,j+1,k
τ
=
(U
ψ
+
U
ψ
+
U
ψ
+
U
ψ
x
y
y
∂t
2h2 x
i,j,k−1 i,j,k−1
+Uzi,j,k ψ i,j,k+1 + U z
ψ
− 6ψ i,j,k ) + η(1 − |ψ i,j,k |2 )ψ i,j,k .
(45
////////// void do_Laplacian(Complex[][][] psi, Complex[][][] Ux, Complex[][][] Uy, Complex[][][] Uz, double dt, double D){ 8
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]); } } } } ////////// 21
境界条件の実装
/// constraints by the applied field and current /// y-z plane ( Bz-Uy-component ) set Ja = 0 when longitudinal for(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 ) longitudinal for(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 ) longitudinal for(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) ); } } 22
今後の課題
•  各種物理量の表示 •  アルゴリズムの改良 •  並列化（高速化） •  ピンニング •  Two-­‐component 系 23