Google ColabのJulia言語で保存系のTDGL方程式(Cahn-Hillard方程式）と非保存系のTDGL方程式(Allen-Cahn方程式）を差分法で計算して動画で比較する。

TDGL方程式(Time Dependent Ginzburg Landau)のオーダーパラメータは保存するものと保存しないものがある。

保存するものはCahn-Hillard方程式とよく呼ばれてこんな形。

∂φ/∂t=∇²（φ³ -φ-∇²φ）

非保存のものはAllen-Cahn方程式と呼ばれることもあってこんな形。

∂φ/∂t=φ-φ³＋∇²φ

ではそれぞれのコードは

保存

using Plots

using Printf

using Random

function main()

    #パラメータ設定

    n = 128

    m = 5000

    ndiv = 20

    dt = 0.01

    dx = 1.0

1

    #初期設定

    X1 = zeros(n + 2, n + 2)

    Y1 = zeros(n + 2, n + 2)

    X2 = zeros(n + 2, n + 2)

    MT = MersenneTwister()

    Random.seed!(MT, 42)

    @inbounds for j in 2:(n + 1)

        @inbounds @simd for i in 2:(n + 1)

            X1[i, j] = 0.1 * (2.0 * rand(MT) - 1.0)

        end

    end

    # 結果を格納する配列

    results = []

    for t in 1:m

        #境界条件

        @inbounds for i in 2:(n + 1)

            X1[1, i] = X1[2, i]

            X1[n + 2, i] = X1[n + 1, i]

            X1[i, 1] = X1[i, 2]

            X1[i, n + 2] = X1[i, n + 1]

        end

        #一つ目のラプラシアン

        @inbounds for j in 2:(n + 1)

            @inbounds @simd for i in 2:(n + 1)

                lapX = (X1[i + 1, j] + X1[i - 1, j] + X1[i, j + 1] + X1[i, j - 1] - 4.0 * X1[i, j]) / (dx * dx)

                Y1[i, j] = X1[i, j] * (X1[i, j]^2 - 1.0) - lapX                                                                              

            end

        end

        #境界条件

        @inbounds for i in 2:(n + 1)    

            Y1[1, i] = Y1[2, i]

            Y1[n + 2, i] = Y1[n + 1, i]

            Y1[i, 1] = Y1[i, 2]

            Y1[i, n + 2] = Y1[i, n + 1]

        end

        #TDGL計算

        @inbounds for j in 2:(n + 1)

            @inbounds @simd for i in 2:(n + 1)

                lapY = (Y1[i + 1, j] + Y1[i - 1, j] + Y1[i, j + 1] + Y1[i, j - 1] - 4.0 * Y1[i, j]) / (dx * dx)

                X2[i, j] = X1[i, j] + dt * lapY                                                                          

            end

        end

        @inbounds for j in 2:(n + 1)

            @inbounds @simd for i in 2:(n + 1)

                X1[i, j] = X2[i, j]

            end

        end

        # ndivステップごとに結果を配列に格納

        if t % ndiv == 0

            push!(results, copy(X1))

        end

    end

    # 計算結果をアニメーションで表示

    anim = @animate for i in 1:length(results)

        heatmap(results[i], title="t = $(@sprintf("%.3f", i * ndiv * dt))", clim=(-1.2, 1.2), size = (950, 800))

    end

    gif(anim, "Cahn-Hillard.gif", fps = 10)    

end

非保存

using Plots

using Printf

using Random

function main()

    #パラメータ設定

    n = 128

    m = 5000

    ndiv = 20

    dt = 0.01

    dx = 1.0

1

    #初期設定

    X1 = zeros(n + 2, n + 2)

    X2 = zeros(n + 2, n + 2)

    MT = MersenneTwister()

    Random.seed!(MT, 42)

    @inbounds for j in 2:(n + 1)

        @inbounds @simd for i in 2:(n + 1)

            X1[i, j] = 0.1 * (2.0 * rand(MT) - 1.0)

        end

    end

    # 結果を格納する配列

    results = []

    for t in 1:m

        #境界条件

        @inbounds for i in 2:(n + 1)

            X1[1, i] = X1[2, i]

            X1[n + 2, i] = X1[n + 1, i]

            X1[i, 1] = X1[i, 2]

            X1[i, n + 2] = X1[i, n + 1]

        end

        #境界条件

        @inbounds for i in 2:(n + 1)    

            X1[1, i] = X1[2, i]

            X1[n + 2, i] = X1[n + 1, i]

            X1[i, 1] = X1[i, 2]

            X1[i, n + 2] = X1[i, n + 1]

        end

        #TDGL計算

        @inbounds for j in 2:(n + 1)

            @inbounds @simd for i in 2:(n + 1)

                lapX = (X1[i + 1, j] + X1[i - 1, j] + X1[i, j + 1] + X1[i, j - 1] - 4.0 * X1[i, j]) / (dx * dx)

                X2[i, j] = X1[i, j] + dt *(X1[i, j] * (1.0 - X1[i, j]^2) + lapX)                                                                          

            end

        end

        @inbounds for j in 2:(n + 1)

            @inbounds @simd for i in 2:(n + 1)

                X1[i, j] = X2[i, j]

            end

        end

        # ndivステップごとに結果を配列に格納

        if t % ndiv == 0

            push!(results, copy(X1))

        end

    end

    # 計算結果をアニメーションで表示

    anim = @animate for i in 1:length(results)

        heatmap(results[i], title="t = $(@sprintf("%.3f", i * ndiv * dt))", clim=(-1.2, 1.2), size = (950, 800))

    end

    gif(anim, "Allen-Cahn.gif", fps = 10)    

end

結果を動画で比較する。

保存

非保存

保存するものは一回構造ができるとほぼ動かないが、非保存のものは島が生成消滅が続く。