wetting/wetting.jl at main · noClaps/wetting · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
using Roots
using DelimitedFiles
using ColorSchemes
import Contour
using Plots
pythonplot()

# Equation whose root defines the homogeneous equilibrium solutions:
function func(beta, mu)
    return (rho) -> rho - (1.0 - rho) * exp(beta * (mu + 5.0 * rho))
end

function problem14(beta_epsilon_wall)
    Lx = 100 # Number of sites along x-axis
    Ly = 40 # Number of sites along y-axis

    beta = 1.2 # Inverse temperature beta*epsilon
    epsilon_wall = beta_epsilon_wall / beta

    mu_coex = -5 / 2
    tol = 1e-12 # Convergence tolerance
    count = 300000 # Upper limit for iterations
    alpha = 0.01 # Mixing parameter

    # Solve equations iteratively:
    conv = 1
    cnt = 1

    minsol = minimum(fzeros(func(beta, mu_coex), 0, 1))
    maxsol = maximum(fzeros(func(beta, mu_coex), 0, 1))

    rho = ones(Lx, Ly) * minsol

    squarewidth = 30
    squareheight = 20
    shape = ones(squarewidth, squareheight) * maxsol
    rho[Int(Lx / 2)-Int(squarewidth / 2)+1:Int(Lx / 2)+Int(squarewidth / 2), 2:squareheight+1] = shape

    rho[:, begin] .= 0

    rho_initial = transpose(rho)
    rho_new = zeros(Lx, Ly)

    while conv >= tol && cnt < count
        cnt += 1
        for i = 1:Lx
            for j = 2:Ly-1
                # Handle the periodic boundaries for x and y:
                left = i == 1 ? Lx : i - 1
                right = i == Lx ? 1 : i + 1
                down = j == 1 ? Ly : j - 1
                up = j == Ly ? 1 : j + 1

                v_j = -epsilon_wall * (j - 1)^(-3)
                rho_new[i, j] = (1 - rho[i, j]) * exp(beta * (rho[i, down] + rho[i, up] + rho[left, j] + rho[right, j] + (1 / 4) * (rho[left, down] + rho[right, down] + rho[left, up] + rho[right, up]) + mu_coex - v_j))
            end
        end

        conv = sum((rho - rho_new) .^ 2) # Compute the convergence parameter.
        rho = alpha * rho_new + (1 - alpha) * rho # Mix the new and old solutions.

        # Normalization step
        N = sum(rho_initial)
        N_current = sum(rho)
        rho *= N / N_current

        rho[:, begin] .= 0
        rho[:, end] .= minsol
    end

    rho = transpose(rho)

    writedlm("problem14-$beta_epsilon_wall.csv", rho, ",")
end

# for beta = 0.1:0.1:2.0
#     problem14(beta)
# end

rho = readdlm("problem14-0.5.csv", ',')
x = 1:size(rho, 1)
y = 1:size(rho, 2)

colorscheme = cgrad(ColorScheme(get(ColorSchemes.jet, range(0.2, 0.7, 256))))
# heatmap(rho, color=colorscheme)
contourf(rho, color=colorscheme, title="Contour plot", xlabel="x positions", ylabel="y positions", colorbar_title="Density gradient", size=(100, 40) .* 12.5, levels=256)

level = Contour.levels(Contour.contours(x, y, rho, 1))[1]
line = first(Contour.lines(level))
ys, xs = Contour.coordinates(line)

indices = []
grads = []

for i = 1:length(xs)-1
    gradient = (ys[i+1] - ys[i]) / (xs[i+1] - xs[i])
    angle = 180 - rad2deg(atan(abs(gradient)))

    if floor(angle) == 115 || ceil(angle) == 180 - 115
        push!(indices, i)
        push!(grads, gradient)
    end
end

for i = eachindex(indices)
    index = indices[i]
    grad = grads[i]

    local x = grad > 0 ? (xs[index]:xs[index]+5) : (xs[index]-5:xs[index])
    local y = (grad .* (x .- xs[index])) .+ ys[index]
    plot!(x, y, color=:black, legend=false)
end

gui()
readline()