時間依存Ginzburg-Landau方程式を解いてみる

Juliaで物理シミュレーションシリーズ。

Juliaで時間発展の方程式を解く練習として、超伝導をシミュレーションするための方程式である時間依存Ginzburg-Landau方程式（通称tdGL、時間依存GL方程式）を解いてみましょう。Juliaのバージョンは1.0.0です。
tdGL方程式を解くことで超伝導体で流せる電流の最大値などが工学分野において計算されていたりします。
https://annex.jsap.or.jp/fluxoid/en/img/FDTD_TDGL_02.pdf

方程式

方程式は

\frac{\partial \Delta({\bf r})}{\partial t} 
= - \left[
\left(-i \nabla + {\bf A}({\bf r}) \right)^2
 \Delta({\bf r}) + (1-T)(|\Delta({\bf r})|^2 -1) \Delta({\bf r}) 
\right]

です。
この方程式を差分化して解いてみます。解き方としては、
http://zvine-ap.eng.hokudai.ac.jp/~asano/pdf2/nishi2007b.pdf
を参考にしました。
ベクトルポテンシャルはMaxwell方程式を解いて求めていません。
ベクトルポテンシャルは

{\bf A}({\bf r})  = \frac{1}{2} {\bf H} \times {\bf r}

で与えられています。なお、磁場は${\bf H} = (0,0,H_z)$です。

コード

コードは以下の通りです。

module TDGL
    using Plots
    struct parameters
        L::Array{Float64,1}
        N::Array{Int64,1}
        dr::Array{Float64,1}               
        Hz::Float64
        δt :: Float64
        Maxstep ::Int64
        filename ::String
    end

    function init_system(L,N,Hz;δt = 0.01,Maxstep=1000,name="GL.dat")
        dr = L ./ (N .-1)
        return parameters(L,N,dr,Hz,δt,Maxstep,name)
    end

    function test()
        L = [25.0,25.0]
        N = [51,51]
        Hz = 0.1
        param = init_system(L,N,Hz,Maxstep = 40000)
        Δ0 = zeros(ComplexF64,N[1],N[2])
        #Δ0 = rand(N[1],N[2])+rand(N[1],N[2])*im
        for ix=1:N[1]
            for iy=1:N[2]
                r = ixiy2r(ix,iy,param) .-5
#                Δ0[ix,iy] = 0.9*exp(2π*im*rand()) #atan(r[2],r[1]))
            end
        end
        Δ0 = 0.9*ones(ComplexF64,N[1],N[2])#.*exp.(2π*im*rand(N[1],N[2]))
#        Δ0[5,5] += im
        T = 0.2#10.0#0.2
        tdGLsimulation(Δ0,param,T)
    end

    Ax(r,param) = -param.Hz*r[2]/2
    Ay(r,param) = param.Hz*r[1]/2

    function ixiy2r(ix,iy,param)
        rx = (ix-1)*param.dr[1]-param.L[1]/2
        ry = (iy-1)*param.dr[2]-param.L[2]/2
        return rx,ry
    end

    function tdGLsimulation(Δ0,param,T)
        Nx = param.N[1]
        Ny = param.N[2]
        δt = param.δt
        name = param.filename
        Δ = copy(Δ0)        
        dΔdt = zeros(ComplexF64,Nx,Ny)
        count = 0
        maps = @animate for i=1:param.Maxstep

            calc_dΔdt!(dΔdt,Δ,param,T)
            Δ += δt*dΔdt            
            heatmap(1:Nx, 1:Ny, abs.(Δ),aspect_ratio=:equal)   
            if i % 40 ==1
                count += 1
                savefig("./files2/git_"*string(count)) 
            end
            println(i,"\t",Δ[div(Nx,2),div(Ny,2)],"\t",abs(Δ[div(Nx,2),div(Ny,2)]),"\t",dΔdt[div(Nx,2),div(Ny,2)])        
            end #every 100
        gif(maps, "./"*name, fps = 15)  
    end

    function calc_dΔdt!(dΔdt,Δ,param,T)
        Nx = param.N[1]
        Ny = param.N[2]
        dr = param.dr
        for ix=1:Nx
            for iy=1:Ny
                dΔdt[ix,iy] = (1-T)*(1-abs(Δ[ix,iy])^2)*Δ[ix,iy]

                for dx =-1:1
                    jx = ix + dx
                    jy = iy
                    if jx > Nx || jx < 1
                        t = 0
                    else
                        if dx == 0
                            t = -2*Δ[ix,iy]/dr[1]^2
                        else
                            r = ixiy2r(jx,jy,param)
                            t = Δ[jx,jy]*exp(dx*im*dr[1]*Ax(r,param))/dr[1]^2
                        end
                    end
                    dΔdt[ix,iy] += t
                end

                for dy =-1:1
                    jx = ix 
                    jy = iy + dy
                    if jy > Ny || jy < 1
                        t = 0
                    else
                        if dy == 0
                            t = -2*Δ[ix,iy]/dr[2]^2
                        else
                            r = ixiy2r(jx,jy,param)
                            t = Δ[jx,jy]*exp(dy*im*dr[2]*Ay(r,param))/dr[2]^2
                        end
                    end
                    dΔdt[ix,iy] += t
                end                
            end
        end        
    end
end

このモジュールにはtest関数をつけているので、

using .TDGL
TDGL.test()

で実行できます。

計算結果

本当はJuliaの@animateを使っているのでそちらでGIFが作成されるはずだったのですが、手元のノートPCにffmpegがうまく入らなかったので、savefigしたpngファイルをあとでフリーソフトでGIFアニメにしました。

その結果が

となります。四角い形状の超伝導体に磁場がかかっているため、量子化された磁場（量子化磁束）が試料内に侵入します。磁場を増やしていくと入る磁束の数が増えていき、ぎゅうぎゅう詰めになってもう入らないような状態になると、超伝導状態が壊れ普通の金属に戻ります。