Introduction
Introduction Statistics Contact Development Disclaimer Help
tbreakup-ocean-atmosphere.jl - seaice-experiments - sea ice experiments using G…
git clone git://src.adamsgaard.dk/seaice-experiments
Log
Files
Refs
README
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 spa…
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=fal…
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, …
151 fixed=true, verbose=false)
152 end
153
154 dx = 2.*r_walls*sin(atan((Lx - Lx_constriction)/(L[2] - Ly_constrict…
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_constr…
159 r_walls, length(x))
160 for i in 1:length(x)
161 Granular.addGrainCylindrical!(sim, [x[i], y[i]], r_walls, h, fix…
162 verbose=false)
163 end
164
165 ## NE diagonal
166 x = (L[1] - r_walls):(-dx):((Lx - Lx_constriction)/2. + Lx_constrict…
167 y = linspace(L[2] - r_walls, (L[2] - Ly_constriction)/2. + Ly_constr…
168 r_walls, length(x))
169 for i in 1:length(x)
170 Granular.addGrainCylindrical!(sim, [x[i], y[i]], r_walls, h, fix…
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, fix…
179 verbose=false)
180 end
181
182 ## SE diagonal
183 x = (L[1] - r_walls):(-dx):((Lx - Lx_constriction)/2. + Lx_constrict…
184 y = linspace(r_walls, (L[2] - Ly_constriction)/2. - r_walls, length(…
185 for i in 1:length(x)
186 Granular.addGrainCylindrical!(sim, [x[i], y[i]], r_walls, h, fix…
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_f…
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.)) …
209
210 elseif sim.ocean.yq[1, j] < Ly_constriction+(L[2] - Ly_constrict…
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)…
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 derivati…
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 f…
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_fl…
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.ocea…
278 r_rand, n_iter=…
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, verbo…
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 durat…
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_st…
343 verbose=verbose, show_file_output=ve…
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 cha…
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…
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" * strin…
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()
You are viewing proxied material from mx1.adamsgaard.dk. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.