FDTD法によるTDGLシミュレーション リンク変数法によるシミュレーション,および オブジェクト指向プログラミングによる込み入った数式の実装 有明高専 電子情報工学科 松野哲也 応用物理学会2015秋(9月13日(日)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 共変微分演算子(2次元) 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 2次元 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 数値シミュレーション • 動画:2次元 横磁界 • 動画:3次元 横磁界,縦磁界,リング • リアルタイムデモ 2次元(40x80),3次元(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
© Copyright 2025 ExpyDoc