simulation.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()