 |
tsimulation.jl - seaice-experiments - sea ice experiments using Granular.jl |
 |
git clone git://src.adamsgaard.dk/seaice-experiments |
 |
Log |
|
Files |
|
Refs |
|
README |
|
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 spa… |
|
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_wa… |
|
150 r_walls, h, fixed=true, verbose=fa… |
|
151 |
|
152 Granular.addGrainCylindrical!(sim, |
|
153 [Lx/2. + Lx_constriction/2. + r_wa… |
|
154 r_walls, h, fixed=true, verbose=fa… |
|
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, fix… |
|
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, fix… |
|
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_f… |
|
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)… |
|
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 derivati… |
|
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 f… |
|
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_fl… |
|
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.ocea… |
|
261 r_rand, n_iter=… |
|
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 durat… |
|
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_st… |
|
339 verbose=verbose, show_file_output=ve… |
|
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 cha… |
|
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_particle… |
|
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_m… |
|
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.… |
|
399 ix, iy, |
|
400 r_rand, n_i… |
|
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_norm… |
|
410 contact_stiffness_tang… |
|
411 k_t, |
|
412 contact_viscosity_norm… |
|
413 gamma_n, |
|
414 contact_viscosity_tang… |
|
415 gamma_t, |
|
416 contact_static_frictio… |
|
417 contact_dynamic_fricti… |
|
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() |
