tjamming.jl - seaice-experiments - sea ice experiments using Granular.jl
(HTM) git clone git://src.adamsgaard.dk/seaice-experiments
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
tjamming.jl (15612B)
---
1 #!/usr/bin/env julia
2 import Granular
3 import ArgParse
4 import Plots
5
6 verbose = false
7
8 function parse_command_line()
9 s = ArgParse.ArgParseSettings()
10 ArgParse.@add_arg_table s begin
11 "--width"
12 help = "strait width [m]"
13 arg_type = Float64
14 default = 5e3
15 "--E"
16 help = "Young's modulus [Pa]"
17 arg_type = Float64
18 default = 0.
19 "--nu"
20 help = "Poisson's ratio [-]"
21 arg_type = Float64
22 default = 0.
23 "--k_n"
24 help = "normal stiffness [N/m]"
25 arg_type = Float64
26 default = 1e7
27 "--k_t"
28 help = "tangential stiffness [N/m]"
29 arg_type = Float64
30 default = 1e7
31 "--gamma_n"
32 help = "normal viscosity [N/(m/s)]"
33 arg_type = Float64
34 default = 0.
35 "--gamma_t"
36 help = "tangential viscosity [N/(m/s)]"
37 arg_type = Float64
38 default = 0.
39 "--mu_s"
40 help = "static friction coefficient [-]"
41 arg_type = Float64
42 default = 0.5
43 "--mu_d"
44 help = "dynamic friction coefficient [-]"
45 arg_type = Float64
46 default = 0.5
47 "--tensile_strength"
48 help = "the maximum tensile strength [Pa]"
49 arg_type = Float64
50 default = 0.
51 "--r_min"
52 help = "minimum grain radius [m]"
53 arg_type = Float64
54 default = 1e3
55 "--r_max"
56 help = "maximum grain radius [m]"
57 arg_type = Float64
58 default = 1e3
59 "--rotating"
60 help = "allow the grains to rotate"
61 arg_type = Bool
62 default = true
63 "--ocean_vel_fac"
64 help = "ocean velocity factor [-]"
65 arg_type = Float64
66 default = 2e4
67 "--total_hours"
68 help = "hours of simulation time [h]"
69 arg_type = Float64
70 default = 6.
71 "--nruns"
72 help = "number of runs in ensemble"
73 arg_type = Int
74 default = 1
75 "id"
76 help = "simulation id"
77 required = true
78 end
79 return ArgParse.parse_args(s)
80 end
81
82 function report_args(parsed_args)
83 println("Parsed args:")
84 for (arg,val) in parsed_args
85 println(" $arg => $val")
86 end
87 end
88
89 function run_simulation(id::String,
90 width::Float64,
91 E::Float64,
92 nu::Float64,
93 k_n::Float64,
94 k_t::Float64,
95 gamma_n::Float64,
96 gamma_t::Float64,
97 mu_s::Float64,
98 mu_d::Float64,
99 tensile_strength::Float64,
100 r_min::Float64,
101 r_max::Float64,
102 rotating::Bool,
103 ocean_vel_fac::Float64,
104 total_hours::Float64,
105 seed::Int)
106
107 info("## EXPERIMENT: " * id * " ##")
108 sim = Granular.createSimulation(id=id)
109
110 Lx = 50.e3
111 Lx_constriction = width
112 L = [Lx, Lx*1.5, 1e3]
113 Ly_constriction = 20e3
114 dx = r_max*2.
115
116 n = [Int(ceil(L[1]/dx/2.)), Int(ceil(L[2]/dx/2.)), 2]
117
118 # Initialize confining walls, which are grains that are fixed in space
119 r = r_max
120 r_min = r/4.
121 h = 1.
122 r_walls = r_min
123
124 ## N-S segments
125 for y in linspace((L[2] - Ly_constriction)/2.,
126 Ly_constriction + (L[2] - Ly_constriction)/2.,
127 Int(round(Ly_constriction/(r_walls*2))))
128 Granular.addGrainCylindrical!(sim, [(Lx - Lx_constriction)/2., y],
129 r_walls, h, fixed=true, verbose=false)
130 end
131 for y in linspace((L[2] - Ly_constriction)/2.,
132 Ly_constriction + (L[2] - Ly_constriction)/2.,
133 Int(round(Ly_constriction/(r_walls*2))))
134 Granular.addGrainCylindrical!(sim, [Lx_constriction + (L[1] -
135 Lx_constriction)/2., y], r_walls, h,
136 fixed=true, verbose=false)
137 end
138
139 dx = 2.*r_walls*sin(atan((Lx - Lx_constriction)/(L[2] - Ly_constriction)))
140
141 ## NW diagonal
142 x = r_walls:dx:((Lx - Lx_constriction)/2.)
143 y = linspace(L[2] - r_walls, (L[2] - Ly_constriction)/2. + Ly_constriction +
144 r_walls, length(x))
145 for i in 1:length(x)
146 Granular.addGrainCylindrical!(sim, [x[i], y[i]], r_walls, h, fixed=true,
147 verbose=false)
148 end
149
150 ## NE diagonal
151 x = (L[1] - r_walls):(-dx):((Lx - Lx_constriction)/2. + Lx_constriction)
152 y = linspace(L[2] - r_walls, (L[2] - Ly_constriction)/2. + Ly_constriction +
153 r_walls, length(x))
154 for i in 1:length(x)
155 Granular.addGrainCylindrical!(sim, [x[i], y[i]], r_walls, h, fixed=true,
156 verbose=false)
157 end
158
159 ## SW diagonal
160 x = r_walls:dx:((Lx - Lx_constriction)/2.)
161 y = linspace(r, (L[2] - Ly_constriction)/2. - r_walls, length(x))
162 for i in 1:length(x)
163 Granular.addGrainCylindrical!(sim, [x[i], y[i]], r_walls, h, fixed=true,
164 verbose=false)
165 end
166
167 ## SE diagonal
168 x = (L[1] - r_walls):(-dx):((Lx - Lx_constriction)/2. + Lx_constriction)
169 y = linspace(r_walls, (L[2] - Ly_constriction)/2. - r_walls, length(x))
170 for i in 1:length(x)
171 Granular.addGrainCylindrical!(sim, [x[i], y[i]], r_walls, h, fixed=true,
172 verbose=false)
173 end
174
175 n_walls = length(sim.grains)
176 info("added $(n_walls) fixed grains as walls")
177
178 # Initialize ocean
179 sim.ocean = Granular.createRegularOceanGrid(n, L, name="poiseuille_flow")
180
181 # Determine stream function value for all grid cells, row by row.
182 # At y=0 and y=Ly: psi = a*ocean_vel_fac*(x - Lx/2.).^3 with a = 1.
183 # The value of a increases when strait is narrow, and psi values are
184 # constant outside domain>
185 psi = similar(sim.ocean.v[:, :, 1, 1])
186 Granular.sortGrainsInGrid!(sim, sim.ocean)
187 for j=1:size(psi, 2)
188
189 # Check width of domain in the current row
190 if sim.ocean.yq[1, j] < (L[2] - Ly_constriction)/2.
191 # lower triangle
192 W = (Lx - width)*
193 (1. - sim.ocean.yq[1, j]/((L[2] - Ly_constriction)/2.)) + width
194
195 elseif sim.ocean.yq[1, j] < Ly_constriction+(L[2] - Ly_constriction)/2.
196 # strait
197 W = width
198
199 else
200 y_min = Ly_constriction + (L[2] - Ly_constriction)/2.
201
202 # upper triangle
203 W = (Lx - width)*(sim.ocean.yq[1, j] - y_min)/(L[2] - y_min) + width
204 end
205
206 # transform [Lx/2 - W/2; Lx/2 + W/2] to [-2;2]
207 x_ = (sim.ocean.xq[:, j] - (Lx/2. - W/2.))/W*4. - 2.
208
209 psi[:, j] = tanh(x_)
210 end
211
212 # determine ocean velocities (u and v) from stream function derivatives
213 for i=1:size(psi, 1)
214 if i == 1
215 sim.ocean.v[i, :, 1, 1] = -(psi[i+1, :] - psi[i, :])./
216 (sim.ocean.xq[i+1, :] - sim.ocean.xq[i, :])
217 elseif i == size(psi, 1)
218 sim.ocean.v[i, :, 1, 1] = -(psi[i, :] - psi[i-1, :])./
219 (sim.ocean.xq[i, :] - sim.ocean.xq[i-1, :])
220 else
221 sim.ocean.v[i, :, 1, 1] = -(psi[i+1, :] - psi[i-1, :])./
222 (sim.ocean.xq[i+1, :] - sim.ocean.xq[i-1, :])
223 end
224 end
225
226 for j=1:size(psi, 2)
227 if j == 1
228 sim.ocean.u[:, j, 1, 1] = (psi[:, j+1] - psi[:, j])./
229 (sim.ocean.yq[:, j+1] - sim.ocean.yq[:, j])
230 elseif j == size(psi, 2)
231 sim.ocean.u[:, j, 1, 1] = (psi[:, j] - psi[:, j-1])./
232 (sim.ocean.yq[:, j] - sim.ocean.yq[:, j-1])
233 else
234 sim.ocean.u[:, j, 1, 1] = (psi[:, j+1] - psi[:, j-1])./
235 (sim.ocean.yq[:, j+1] - sim.ocean.yq[:, j-1])
236 end
237 end
238 sim.ocean.h[:,:,1,1] = psi # this field is unused; use for stream function
239 sim.ocean.u *= ocean_vel_fac
240 sim.ocean.v *= ocean_vel_fac
241
242 # Constant velocities along y:
243 #sim.ocean.v[:, :, 1, 1] = ocean_vel_fac*((sim.ocean.xq - Lx/2.).^2 -
244 # Lx^2./4.)
245
246
247 # Initialize grains in wedge north of the constriction
248 iy = 1
249 dy = sqrt((2.*r_walls)^2. - dx^2.)
250 spacing_to_boundaries = 4.*r
251 floe_padding = .5*r
252 noise_amplitude = floe_padding
253 Base.Random.srand(seed)
254 for y in (L[2] - r - noise_amplitude):(-(2.*r + floe_padding)):((L[2] -
255 Ly_constriction)/2. + Ly_constriction)
256 for x in (r + noise_amplitude):(2.*r + floe_padding):(Lx - r -
257 noise_amplitude)
258 if iy % 2 == 0
259 x += 1.5*r
260 end
261
262 x_ = x + noise_amplitude*(0.5 - Base.Random.rand())
263 y_ = y + noise_amplitude*(0.5 - Base.Random.rand())
264
265 if y_ < -dy/dx*x_ + L[2] + spacing_to_boundaries
266 continue
267 end
268
269 if y_ < dy/dx*x_ + (L[2] - dy/dx*Lx) + spacing_to_boundaries
270 continue
271 end
272
273 r_rand = r_min + Base.Random.rand()*(r - r_min)
274 Granular.addGrainCylindrical!(sim, [x_, y_], r_rand, h,
275 verbose=false)
276 end
277 iy += 1
278 end
279 n = length(sim.grains) - n_walls
280 info("added $(n) grains")
281
282 # Remove old simulation files
283 Granular.removeSimulationFiles(sim)
284
285 for i=1:length(sim.grains)
286 sim.grains[i].youngs_modulus = E
287 sim.grains[i].poissons_ratio = nu
288 sim.grains[i].contact_stiffness_normal = k_n
289 sim.grains[i].contact_stiffness_tangential = k_t
290 sim.grains[i].contact_viscosity_normal = gamma_n
291 sim.grains[i].contact_viscosity_tangential = gamma_t
292 sim.grains[i].contact_static_friction = mu_s
293 sim.grains[i].contact_dynamic_friction = mu_d
294 sim.grains[i].tensile_strength = tensile_strength
295 sim.grains[i].rotating = rotating
296 end
297
298 # Set temporal parameters
299 Granular.setTotalTime!(sim, total_hours*60.*60.)
300 Granular.setOutputFileInterval!(sim, 60.*5.) # output every 5 mins
301 Granular.setTimeStep!(sim, verbose=verbose)
302 Granular.writeVTK(sim, verbose=verbose)
303
304 profile = false
305
306 ice_flux = Float64[]
307 jammed = false
308 it_before_eval = 10
309 time_jammed = 0.
310 time_jammed_criterion = 60.*60. # define as jammed after this duration
311
312 while sim.time < sim.time_total
313
314 # run simulation for it_before_eval time steps
315 for it=1:it_before_eval
316 if sim.time >= sim.time_total*.75 && profile
317 @profile Granular.run!(sim, status_interval=1, single_step=true,
318 verbose=verbose, show_file_output=verbose)
319 Profile.print()
320 profile = false
321 else
322 Granular.run!(sim, status_interval=1, single_step=true,
323 verbose=verbose, show_file_output=verbose)
324 end
325 end
326
327 # assert if the system is jammed by looking at ice-floe mass change in
328 # the number of jammed grains
329 if jammed
330 time_jammed += sim.time_dt*float(it_before_eval)
331 if time_jammed >= 60.*60. # 1 h
332 info("$t s: system jammed for more than " *
333 "$time_jammed_criterion s, stopping simulation")
334 exit()
335 end
336 end
337
338 ice_mass_outside_domain = 0.
339 for icefloe in sim.grains
340 if !icefloe.enabled
341 ice_mass_outside_domain += icefloe.mass
342 end
343 end
344 append!(ice_flux, ice_mass_outside_domain)
345
346 # add new grains from the top
347 for i=1:size(sim.ocean.xh, 1)
348 if sim.ocean.grain_list[i, end] == []
349
350 x, y = Granular.getCellCenterCoordinates(sim.ocean, i,
351 size(sim.ocean.xh, 2))
352
353 x_ = x + noise_amplitude*(0.5 - Base.Random.rand())
354 y_ = y + noise_amplitude*(0.5 - Base.Random.rand())
355 r_rand = r_min + Base.Random.rand()*(r - r_min)
356
357 Granular.addGrainCylindrical!(sim, [x_, y_], r_rand, h,
358 verbose=false,
359 youngs_modulus=E,
360 poissons_ratio=nu,
361 contact_stiffness_normal=k_n,
362 contact_stiffness_tangential=k_t,
363 contact_viscosity_normal=gamma_n,
364 contact_viscosity_tangential=gamma_t,
365 contact_static_friction = mu_s,
366 contact_dynamic_friction = mu_d,
367 tensile_strength = tensile_strength,
368 rotating=rotating)
369 end
370 end
371 end
372
373 t = linspace(0., sim.time_total, length(ice_flux))
374 writedlm(sim.id * "-ice-flux.txt", [t, ice_flux])
375 return t, ice_flux
376 end
377
378 function main()
379 parsed_args = parse_command_line()
380 report_args(parsed_args)
381
382 nruns = parsed_args["nruns"]
383 t = Float64[]
384 t_jam = ones(nruns)*parsed_args["total_hours"]*60.*60.*2.
385 Plots.gr()
386
387 for i=1:nruns
388 seed = i
389 t, ice_flux = run_simulation(parsed_args["id"] * "-seed" * string(seed),
390 parsed_args["width"],
391 parsed_args["E"],
392 parsed_args["nu"],
393 parsed_args["k_n"],
394 parsed_args["k_t"],
395 parsed_args["gamma_n"],
396 parsed_args["gamma_t"],
397 parsed_args["mu_s"],
398 parsed_args["mu_d"],
399 parsed_args["tensile_strength"],
400 parsed_args["r_min"],
401 parsed_args["r_max"],
402 parsed_args["rotating"],
403 parsed_args["ocean_vel_fac"],
404 parsed_args["total_hours"],
405 i)
406 Plots.plot!(t/(60.*60.), ice_flux)
407
408 time_elapsed_while_jammed = 0.
409 for it=(length(t) - 1):-1:1
410 if ice_flux[it] ≈ ice_flux[end]
411 time_elapsed_while_jammed = t[end] - t[it]
412 end
413 end
414
415 if time_elapsed_while_jammed > 60.*60.
416 t_jam[i] = t[end] - time_elapsed_while_jammed
417 info("simulation $i jammed at t = $(t_jam[i]/(60.*60.)) h")
418 end
419
420 end
421 Plots.title!(parsed_args["id"])
422 Plots.xaxis!("Time [h]")
423 Plots.yaxis!("Cumulative ice throughput [kg]")
424
425 Plots.savefig(parsed_args["id"])
426 Plots.closeall()
427
428 jam_fraction = zeros(length(t))
429 for it=1:length(t)
430 for i=1:nruns
431 if t_jam[i] <= t[it]
432 jam_fraction[it] += 1./float(nruns)
433 end
434 end
435 end
436 Plots.plot(t/(60.*60.), jam_fraction)
437 Plots.title!(parsed_args["id"] * ", N = " * string(nruns))
438 Plots.xaxis!("Time [h]")
439 Plots.yaxis!("Fraction of runs jammed [-]")
440 Plots.savefig(parsed_args["id"] * "-jam_fraction.pdf")
441 end
442
443 main()