 |
tjamming.jl - seaice-experiments - sea ice experiments using Granular.jl |
 |
git clone git://src.adamsgaard.dk/seaice-experiments |
 |
Log |
|
Files |
|
Refs |
|
README |
|
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 spa… |
|
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=fal… |
|
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, … |
|
136 fixed=true, verbose=false) |
|
137 end |
|
138 |
|
139 dx = 2.*r_walls*sin(atan((Lx - Lx_constriction)/(L[2] - Ly_constrict… |
|
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_constr… |
|
144 r_walls, length(x)) |
|
145 for i in 1:length(x) |
|
146 Granular.addGrainCylindrical!(sim, [x[i], y[i]], r_walls, h, fix… |
|
147 verbose=false) |
|
148 end |
|
149 |
|
150 ## NE diagonal |
|
151 x = (L[1] - r_walls):(-dx):((Lx - Lx_constriction)/2. + Lx_constrict… |
|
152 y = linspace(L[2] - r_walls, (L[2] - Ly_constriction)/2. + Ly_constr… |
|
153 r_walls, length(x)) |
|
154 for i in 1:length(x) |
|
155 Granular.addGrainCylindrical!(sim, [x[i], y[i]], r_walls, h, fix… |
|
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, fix… |
|
164 verbose=false) |
|
165 end |
|
166 |
|
167 ## SE diagonal |
|
168 x = (L[1] - r_walls):(-dx):((Lx - Lx_constriction)/2. + Lx_constrict… |
|
169 y = linspace(r_walls, (L[2] - Ly_constriction)/2. - r_walls, length(… |
|
170 for i in 1:length(x) |
|
171 Granular.addGrainCylindrical!(sim, [x[i], y[i]], r_walls, h, fix… |
|
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_f… |
|
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.)) … |
|
194 |
|
195 elseif sim.ocean.yq[1, j] < Ly_constriction+(L[2] - Ly_constrict… |
|
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)… |
|
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 derivati… |
|
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 f… |
|
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_ampl… |
|
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 durat… |
|
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_st… |
|
318 verbose=verbose, show_file_output=ve… |
|
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 cha… |
|
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… |
|
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" * strin… |
|
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() |
