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