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