-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathsimulate.jl
170 lines (168 loc) · 8.47 KB
/
simulate.jl
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
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
@testset "Simulations" begin
T = Float64
@testset "Flat Interface" begin
sys = Swalbe.SysConst(Lx=25, Ly=25, param=Swalbe.Taumucs(Tmax=200, tdump=100))
h = Swalbe.run_flat(sys, "CPU", verbos=false)
@test all(h .== 1.0)
@test sum(h) == 25*25
h = Swalbe.run_flat(sys, "CPU", verbos=true)
# One dim version
sys = Swalbe.SysConst_1D(L=25, param=Swalbe.Taumucs(Tmax=200, tdump=100))
h = Swalbe.run_flat(sys, verbos=false)
@test all(h .== 1.0)
@test sum(h) == 25
h = Swalbe.run_flat(sys, verbos=true)
end
@testset "Random Interface" begin
sys = Swalbe.SysConst(Lx=25, Ly=25, param=Swalbe.Taumucs(Tmax=10000, tdump=5000))
h = Swalbe.run_random(sys, "CPU", ϵ=0.1, verbos=false)
difference = maximum(h) - minimum(h)
# Interface should flatten with time, thus difference must be less than the initial ϵ.
@test difference < 0.1
sys = Swalbe.SysConst(Lx=25, Ly=25, param=Swalbe.Taumucs(Tmax=100, tdump=500))
h = Swalbe.run_random(sys, "CPU", ϵ=0.1, verbos=true)
# One dim version
sys = Swalbe.SysConst_1D(L=25, param=Swalbe.Taumucs(Tmax=10000, tdump=5000))
h = Swalbe.run_random(sys, ϵ=0.1, verbos=false)
difference = maximum(h) - minimum(h)
@test difference < 0.02
sys = Swalbe.SysConst_1D(L=25, param=Swalbe.Taumucs(Tmax=100, tdump=500))
h = Swalbe.run_random(sys, ϵ=0.1, verbos=true)
end
@testset "Rayleigh Taylor" begin
sys = Swalbe.SysConst(Lx=100, Ly=100, param=Swalbe.Taumucs(Tmax=1000, tdump=500, g=-0.002))
_, diff = Swalbe.run_rayleightaylor(sys, "CPU", kx=4, ky=5, ϵ=0.01, verbos=true)
@test diff[1] < diff[end]
# TODO: Write a test to check the growthrate
# One dim case
sys = Swalbe.SysConst_1D(L = 512, param=Swalbe.Taumucs(Tmax=1000, tdump=500, g=-0.002))
_, diff = Swalbe.run_rayleightaylor(sys, k=4, ϵ=0.01, verbos=true)
@test diff[1] < diff[end]
end
@testset "Relaxing droplet" begin
rads = 35
sys = Swalbe.SysConst(Lx=150, Ly=150, param=Swalbe.Taumucs(Tmax=10000, δ=3.0))
h, A = Swalbe.run_dropletrelax(sys, "CPU", radius=rads, verbos=false)
# Initial droplet volume
vol = π/3 * rads^3 * (2+cospi(1/6)) * (1-cospi(1/6))^2
R1 = cbrt((rads^3*(2+cospi(1/6))*(1-cospi(1/6))^2)/((2+cospi(1/9))*(1-cospi(1/9))^2))
r1 = sinpi(1/9)*R1
# Identify the droplet, below this height is precursor layer
drop1d = findall(h[75,:] .> 0.055)
# Get the radius and height
droprad = length(drop1d)/2
droph = maximum(h)
# Compute the droplets shape after 10k time steps
vnum = 1/6*π*droph*(3*droprad^2 + droph^2)
# Test that the droplet volume has not changed too much
@test vol ≈ vnum atol = vol/100*10
@test r1 ≈ droprad atol = r1/100*10
@test A[1] < A[end]
sys = Swalbe.SysConst(Lx=150, Ly=150, param=Swalbe.Taumucs(Tmax=100, δ=3.0, tdump=50))
h, A = Swalbe.run_dropletrelax(sys, "CPU", radius=rads, verbos=true)
# One dim case
thresh = 0.05026
sys = Swalbe.SysConst_1D(L = 150, param=Swalbe.Taumucs(Tmax=10000, δ=3.0))
h, A = Swalbe.run_dropletrelax(sys, radius=rads, verbos=false)
# Initial droplet volume
vol = π/3 * rads^3 * (2+cospi(1/6)) * (1-cospi(1/6))^2
R1 = cbrt((rads^3*(2+cospi(1/6))*(1-cospi(1/6))^2)/((2+cospi(1/9))*(1-cospi(1/9))^2))
r1 = sinpi(1/9)*R1
# Identify the droplet, below this height is precursor layer
drop1d = findall(h .> 0.055)
# Drop area (summation of height)
drop_ar = zeros(sys.L)
for i in 1:sys.L
if h[i] > thresh
drop_ar[i] = h[i]
end
end
# Get the radius and height
droprad = length(drop1d)/2
droph = maximum(h)
# Compute the droplets shape after 10k time steps
vnum = 1/6*π*droph*(3*droprad^2 + droph^2)
sph_r = droph/2 + length(findall(h .> thresh))^2/(8*droph)
ar = 2/3 * length(drop1d) * droph # sph_r * π/6 # It's one dimensional so it is the arc length .. sph_r^2/2*(π/9 - sinpi(1/9))
# Test that the droplet volume has not changed too much
@test vol ≈ vnum atol = vol/100*10
# TODO fix these two tests
@test ar ≈ sum(drop_ar) atol = ar/100*10
@test r1 ≈ droprad atol = r1/100*10
@test A[1] < A[end]
sys = Swalbe.SysConst_1D(L = 200, param=Swalbe.Taumucs(Tmax=100, δ=3.0, tdump=50))
h, A = Swalbe.run_dropletrelax(sys, radius=35, verbos=true)
end
@testset "Patterned droplet" begin
sys = Swalbe.SysConst(Lx=150, Ly=150, param=Swalbe.Taumucs(Tmax=10000, δ=3.0))
h = Swalbe.run_dropletpatterned(sys, "CPU", radius=35, θₛ = fill(1/9,sys.Lx, sys.Ly), verbos=false)
# Initial droplet volume
vol = π/3 * 35^3 * (2+cospi(1/6)) * (1-cospi(1/6))^2
R1 = cbrt((35^3*(2+cospi(1/6))*(1-cospi(1/6))^2)/((2+cospi(1/9))*(1-cospi(1/9))^2))
r1 = sinpi(1/9)*R1
# Identify the droplet, below this height is precursor layer
drop1d = findall(h[75,:] .> 0.055)
# Get the radius and height
droprad = length(drop1d)/2
droph = maximum(h)
# Compute the droplets shape after 10k time steps
vnum = 1/6*π*droph*(3*droprad^2 + droph^2)
# Test that the droplet volume has not changed too much
@test vol ≈ vnum atol = vol/100*10
@test r1 ≈ droprad atol = r1/100*10
sys = Swalbe.SysConst(Lx=150, Ly=150, param=Swalbe.Taumucs(Tmax=100, tdump=50, δ=3.0))
h = Swalbe.run_dropletpatterned(sys, "CPU", radius=35, θₛ = fill(1/9,sys.Lx, sys.Ly))
# One dim case
thresh = 0.05026
sys = Swalbe.SysConst_1D(L=150, param=Swalbe.Taumucs(Tmax=10000, δ=3.0))
h = Swalbe.run_dropletpatterned(sys, radius=35, θₛ = fill(1/9,sys.L), verbos=false)
# Initial droplet volume
vol = π/3 * 35^3 * (2+cospi(1/6)) * (1-cospi(1/6))^2
R1 = cbrt((35^3*(2+cospi(1/6))*(1-cospi(1/6))^2)/((2+cospi(1/9))*(1-cospi(1/9))^2))
r1 = sinpi(1/9)*R1
# Identify the droplet, below this height is precursor layer
drop1d = findall(h .> 0.055)
# Drop area (summation of height)
drop_ar = zeros(sys.L)
for i in 1:sys.L
if h[i] > thresh
drop_ar[i] = h[i]
end
end
# Get the radius and height
droprad = length(drop1d)/2
droph = maximum(h)
# Compute the droplets shape after 10k time steps
vnum = 1/6*π*droph*(3*droprad^2 + droph^2)
sph_r = droph/2 + length(findall(h .> thresh))^2/(8*droph)
ar = sph_r^2/2*(π/9 - sinpi(1/9))
# TODO fix me
# @test ar ≈ sum(drop_ar) atol = ar/100*10
@test r1 ≈ droprad atol = r1/100*10
sys = Swalbe.SysConst_1D(L=150, param=Swalbe.Taumucs(Tmax=100, tdump=50, δ=3.0))
h = Swalbe.run_dropletpatterned(sys, radius=35, θₛ = fill(1/9,sys.L))
end
@testset "Sliding droplet" begin
sys = Swalbe.SysConst(Lx=150, Ly=150, param=Swalbe.Taumucs(Tmax=5000, δ=2.0))
# This will return the hydrodynamic moments: h, u, v
mom = Swalbe.run_dropletforced(sys, "CPU", radius=35, fx=1e-4, verbos=false)
# Now check that the center has moved
@test findmax(mom[1])[2][1] ≠ 75
@test findmax(mom[1])[2][2] == 75
# Further test that the velocity is not too large
@test all(mom[2] .< 0.1)
@test all(mom[3] .< 0.1)
sys = Swalbe.SysConst(Lx=150, Ly=150, param=Swalbe.Taumucs(Tmax=100, tdump=50, δ=3.0))
h = Swalbe.run_dropletforced(sys, "CPU", radius=35, fx=1e-4)
# One dim case
sys = Swalbe.SysConst_1D(L=150, param=Swalbe.Taumucs(Tmax=5000, δ=2.0))
# This will return the hydrodynamic moments: h, u, v
mom = Swalbe.run_dropletforced(sys, radius=35, f=1e-4, verbos=false)
# Now check that the center has moved
@test findmax(mom[1])[2] ≠ Int(sys.L/2)
# Further test that the velocity is not too large
@test all(mom[2] .< 0.1)
sys = Swalbe.SysConst_1D(L=150, param=Swalbe.Taumucs(Tmax=100, tdump=50, δ=3.0))
h = Swalbe.run_dropletforced(sys, radius=35, f=1e-4)
end
end