jamming.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 space 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=false) 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, h, 136 fixed=true, verbose=false) 137 end 138 139 dx = 2.*r_walls*sin(atan((Lx - Lx_constriction)/(L[2] - Ly_constriction))) 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_constriction + 144 r_walls, length(x)) 145 for i in 1:length(x) 146 Granular.addGrainCylindrical!(sim, [x[i], y[i]], r_walls, h, fixed=true, 147 verbose=false) 148 end 149 150 ## NE diagonal 151 x = (L[1] - r_walls):(-dx):((Lx - Lx_constriction)/2. + Lx_constriction) 152 y = linspace(L[2] - r_walls, (L[2] - Ly_constriction)/2. + Ly_constriction + 153 r_walls, length(x)) 154 for i in 1:length(x) 155 Granular.addGrainCylindrical!(sim, [x[i], y[i]], r_walls, h, fixed=true, 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, fixed=true, 164 verbose=false) 165 end 166 167 ## SE diagonal 168 x = (L[1] - r_walls):(-dx):((Lx - Lx_constriction)/2. + Lx_constriction) 169 y = linspace(r_walls, (L[2] - Ly_constriction)/2. - r_walls, length(x)) 170 for i in 1:length(x) 171 Granular.addGrainCylindrical!(sim, [x[i], y[i]], r_walls, h, fixed=true, 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_flow") 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.)) + width 194 195 elseif sim.ocean.yq[1, j] < Ly_constriction+(L[2] - Ly_constriction)/2. 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) + width 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 derivatives 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 function 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_amplitude) 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 duration 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_step=true, 318 verbose=verbose, show_file_output=verbose) 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 change in 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, 2)) 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" * string(seed), 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()