breakup.jl (15190B)
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" 52 help = "grain radius [m]" 53 arg_type = Float64 54 default = 1e3 55 "--rotating" 56 help = "allow the grains to rotate" 57 arg_type = Bool 58 default = true 59 "--ocean_vel_fac" 60 help = "ocean velocity factor [-]" 61 arg_type = Float64 62 default = 2e4 63 "--total_hours" 64 help = "hours of simulation time [h]" 65 arg_type = Float64 66 default = 6. 67 "--nruns" 68 help = "number of runs in ensemble" 69 arg_type = Int 70 default = 1 71 "id" 72 help = "simulation id" 73 required = true 74 end 75 return ArgParse.parse_args(s) 76 end 77 78 function report_args(parsed_args) 79 println("Parsed args:") 80 for (arg,val) in parsed_args 81 println(" $arg => $val") 82 end 83 end 84 85 function run_simulation(id::String, 86 width::Float64, 87 E::Float64, 88 nu::Float64, 89 k_n::Float64, 90 k_t::Float64, 91 gamma_n::Float64, 92 gamma_t::Float64, 93 mu_s::Float64, 94 mu_d::Float64, 95 tensile_strength::Float64, 96 r::Float64, 97 rotating::Bool, 98 ocean_vel_fac::Float64, 99 total_hours::Float64, 100 seed::Int) 101 102 info("## EXPERIMENT: " * id * " ##") 103 sim = Granular.createSimulation(id=id) 104 105 Lx = 50.e3 106 Lx_constriction = width 107 L = [Lx, Lx*1.5, 1e3] 108 Ly_constriction = 20e3 109 dx = r*2. 110 111 n = [Int(ceil(L[1]/dx/2.)), Int(ceil(L[2]/dx/2.)), 2] 112 113 # Initialize confining walls, which are grains that are fixed in space 114 h = 1. 115 r_walls = r/2. 116 117 ## N-S segments 118 for y in linspace((L[2] - Ly_constriction)/2., 119 Ly_constriction + (L[2] - Ly_constriction)/2., 120 Int(round(Ly_constriction/(r_walls*2)))) 121 Granular.addGrainCylindrical!(sim, [(Lx - Lx_constriction)/2., y], 122 r_walls, h, fixed=true, verbose=false) 123 end 124 for y in linspace((L[2] - Ly_constriction)/2., 125 Ly_constriction + (L[2] - Ly_constriction)/2., 126 Int(round(Ly_constriction/(r_walls*2)))) 127 Granular.addGrainCylindrical!(sim, [Lx_constriction + (L[1] - 128 Lx_constriction)/2., y], r_walls, h, 129 fixed=true, verbose=false) 130 end 131 132 dx = 2.*r_walls*sin(atan((Lx - Lx_constriction)/(L[2] - Ly_constriction))) 133 134 ## NW diagonal 135 x = r_walls:dx:((Lx - Lx_constriction)/2.) 136 y = linspace(L[2] - r_walls, (L[2] - Ly_constriction)/2. + Ly_constriction + 137 r_walls, length(x)) 138 for i in 1:length(x) 139 Granular.addGrainCylindrical!(sim, [x[i], y[i]], r_walls, h, fixed=true, 140 verbose=false) 141 end 142 143 ## NE diagonal 144 x = (L[1] - r_walls):(-dx):((Lx - Lx_constriction)/2. + Lx_constriction) 145 y = linspace(L[2] - r_walls, (L[2] - Ly_constriction)/2. + Ly_constriction + 146 r_walls, length(x)) 147 for i in 1:length(x) 148 Granular.addGrainCylindrical!(sim, [x[i], y[i]], r_walls, h, fixed=true, 149 verbose=false) 150 end 151 152 ## SW diagonal 153 x = r_walls:dx:((Lx - Lx_constriction)/2.) 154 y = linspace(r, (L[2] - Ly_constriction)/2. - r_walls, length(x)) 155 for i in 1:length(x) 156 Granular.addGrainCylindrical!(sim, [x[i], y[i]], r_walls, h, fixed=true, 157 verbose=false) 158 end 159 160 ## SE diagonal 161 x = (L[1] - r_walls):(-dx):((Lx - Lx_constriction)/2. + Lx_constriction) 162 y = linspace(r_walls, (L[2] - Ly_constriction)/2. - 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 n_walls = length(sim.grains) 169 info("added $(n_walls) fixed grains as walls") 170 171 # Initialize ocean 172 sim.ocean = Granular.createRegularOceanGrid(n, L, name="poiseuille_flow") 173 174 # Determine stream function value for all grid cells, row by row. 175 # At y=0 and y=Ly: psi = a*ocean_vel_fac*(x - Lx/2.).^3 with a = 1. 176 # The value of a increases when strait is narrow, and psi values are 177 # constant outside domain> 178 psi = similar(sim.ocean.v[:, :, 1, 1]) 179 Granular.sortGrainsInGrid!(sim, sim.ocean) 180 for j=1:size(psi, 2) 181 182 # Check width of domain in the current row 183 if sim.ocean.yq[1, j] < (L[2] - Ly_constriction)/2. 184 # lower triangle 185 W = (Lx - width)* 186 (1. - sim.ocean.yq[1, j]/((L[2] - Ly_constriction)/2.)) + width 187 188 elseif sim.ocean.yq[1, j] < Ly_constriction+(L[2] - Ly_constriction)/2. 189 # strait 190 W = width 191 192 else 193 y_min = Ly_constriction + (L[2] - Ly_constriction)/2. 194 195 # upper triangle 196 W = (Lx - width)*(sim.ocean.yq[1, j] - y_min)/(L[2] - y_min) + width 197 end 198 199 # transform [Lx/2 - W/2; Lx/2 + W/2] to [-2;2] 200 x_ = (sim.ocean.xq[:, j] - (Lx/2. - W/2.))/W*4. - 2. 201 202 psi[:, j] = tanh(x_) 203 end 204 205 # determine ocean velocities (u and v) from stream function derivatives 206 for i=1:size(psi, 1) 207 if i == 1 208 sim.ocean.v[i, :, 1, 1] = -(psi[i+1, :] - psi[i, :])./ 209 (sim.ocean.xq[i+1, :] - sim.ocean.xq[i, :]) 210 elseif i == size(psi, 1) 211 sim.ocean.v[i, :, 1, 1] = -(psi[i, :] - psi[i-1, :])./ 212 (sim.ocean.xq[i, :] - sim.ocean.xq[i-1, :]) 213 else 214 sim.ocean.v[i, :, 1, 1] = -(psi[i+1, :] - psi[i-1, :])./ 215 (sim.ocean.xq[i+1, :] - sim.ocean.xq[i-1, :]) 216 end 217 end 218 219 for j=1:size(psi, 2) 220 if j == 1 221 sim.ocean.u[:, j, 1, 1] = (psi[:, j+1] - psi[:, j])./ 222 (sim.ocean.yq[:, j+1] - sim.ocean.yq[:, j]) 223 elseif j == size(psi, 2) 224 sim.ocean.u[:, j, 1, 1] = (psi[:, j] - psi[:, j-1])./ 225 (sim.ocean.yq[:, j] - sim.ocean.yq[:, j-1]) 226 else 227 sim.ocean.u[:, j, 1, 1] = (psi[:, j+1] - psi[:, j-1])./ 228 (sim.ocean.yq[:, j+1] - sim.ocean.yq[:, j-1]) 229 end 230 end 231 sim.ocean.h[:,:,1,1] = psi # this field is unused; use for stream function 232 sim.ocean.u *= ocean_vel_fac 233 sim.ocean.v *= ocean_vel_fac 234 235 # Constant velocities along y: 236 #sim.ocean.v[:, :, 1, 1] = ocean_vel_fac*((sim.ocean.xq - Lx/2.).^2 - 237 # Lx^2./4.) 238 239 240 # Initialize grains in wedge north of the constriction 241 iy = 1 242 dy = sqrt((2.*r_walls)^2. - dx^2.) 243 spacing_to_boundaries = 2.*r 244 #floe_padding = .5*r 245 floe_padding = 0. 246 noise_amplitude = floe_padding 247 Base.Random.srand(seed) 248 y_spacing = 2.*r*cos(pi/6.) + floe_padding 249 for y in (L[2] - r - noise_amplitude):(-y_spacing):((L[2] - 250 Ly_constriction)/2. + Ly_constriction) 251 for x in (r + noise_amplitude):(2.*r + floe_padding):(Lx - r - 252 noise_amplitude) 253 if iy % 2 == 0 254 x += r 255 end 256 257 x_ = x + noise_amplitude*(0.5 - Base.Random.rand()) 258 y_ = y + noise_amplitude*(0.5 - Base.Random.rand()) 259 260 if y_ < -dy/dx*x_ + L[2] + spacing_to_boundaries 261 continue 262 end 263 264 if y_ < dy/dx*x_ + (L[2] - dy/dx*Lx) + spacing_to_boundaries 265 continue 266 end 267 268 Granular.addGrainCylindrical!(sim, [x_, y_], r, h, verbose=false) 269 end 270 iy += 1 271 end 272 n = length(sim.grains) - n_walls 273 info("added $(n) grains") 274 275 # Remove old simulation files 276 Granular.removeSimulationFiles(sim) 277 278 for i=1:length(sim.grains) 279 sim.grains[i].youngs_modulus = E 280 sim.grains[i].poissons_ratio = nu 281 sim.grains[i].contact_stiffness_normal = k_n 282 sim.grains[i].contact_stiffness_tangential = k_t 283 sim.grains[i].contact_viscosity_normal = gamma_n 284 sim.grains[i].contact_viscosity_tangential = gamma_t 285 sim.grains[i].contact_static_friction = mu_s 286 sim.grains[i].contact_dynamic_friction = mu_d 287 sim.grains[i].tensile_strength = tensile_strength 288 sim.grains[i].rotating = rotating 289 end 290 291 # Set temporal parameters 292 Granular.setTotalTime!(sim, total_hours*60.*60.) 293 Granular.setOutputFileInterval!(sim, 60.*5.) # output every 5 mins 294 Granular.setTimeStep!(sim, verbose=verbose) 295 Granular.writeVTK(sim, verbose=verbose) 296 297 profile = false 298 299 ice_flux = Float64[] 300 jammed = false 301 it_before_eval = 10 302 time_jammed = 0. 303 time_jammed_criterion = 60.*60. # define as jammed after this duration 304 305 while sim.time < sim.time_total 306 307 # run simulation for it_before_eval time steps 308 for it=1:it_before_eval 309 if sim.time >= sim.time_total*.75 && profile 310 @profile Granular.run!(sim, status_interval=1, single_step=true, 311 verbose=verbose, show_file_output=verbose) 312 Profile.print() 313 profile = false 314 else 315 Granular.run!(sim, status_interval=1, single_step=true, 316 verbose=verbose, show_file_output=verbose) 317 end 318 end 319 320 # assert if the system is jammed by looking at ice-floe mass change in 321 # the number of jammed grains 322 if jammed 323 time_jammed += sim.time_dt*float(it_before_eval) 324 if time_jammed >= 60.*60. # 1 h 325 info("$t s: system jammed for more than " * 326 "$time_jammed_criterion s, stopping simulation") 327 exit() 328 end 329 end 330 331 ice_mass_outside_domain = 0. 332 for icefloe in sim.grains 333 if !icefloe.enabled 334 ice_mass_outside_domain += icefloe.mass 335 end 336 end 337 append!(ice_flux, ice_mass_outside_domain) 338 339 # add new grains from the top 340 for i=1:size(sim.ocean.xh, 1) 341 if sim.ocean.grain_list[i, end] == [] 342 343 x, y = Granular.getCellCenterCoordinates(sim.ocean, i, 344 size(sim.ocean.xh, 2)) 345 346 x_ = x + noise_amplitude*(0.5 - Base.Random.rand()) 347 y_ = y + noise_amplitude*(0.5 - Base.Random.rand()) 348 349 Granular.addGrainCylindrical!(sim, [x_, y_], r, h, verbose=false, 350 youngs_modulus=E, 351 poissons_ratio=nu, 352 contact_stiffness_normal=k_n, 353 contact_stiffness_tangential=k_t, 354 contact_viscosity_normal=gamma_n, 355 contact_viscosity_tangential=gamma_t, 356 contact_static_friction = mu_s, 357 contact_dynamic_friction = mu_d, 358 tensile_strength = tensile_strength, 359 rotating=rotating) 360 end 361 end 362 end 363 364 t = linspace(0., sim.time_total, length(ice_flux)) 365 writedlm(sim.id * "-ice-flux.txt", [t, ice_flux]) 366 return t, ice_flux 367 end 368 369 function main() 370 parsed_args = parse_command_line() 371 report_args(parsed_args) 372 373 nruns = parsed_args["nruns"] 374 t = Float64[] 375 t_jam = ones(nruns)*parsed_args["total_hours"]*60.*60.*2. 376 Plots.gr() 377 378 for i=1:nruns 379 seed = i 380 t, ice_flux = run_simulation(parsed_args["id"] * "-seed" * string(seed), 381 parsed_args["width"], 382 parsed_args["E"], 383 parsed_args["nu"], 384 parsed_args["k_n"], 385 parsed_args["k_t"], 386 parsed_args["gamma_n"], 387 parsed_args["gamma_t"], 388 parsed_args["mu_s"], 389 parsed_args["mu_d"], 390 parsed_args["tensile_strength"], 391 parsed_args["r"], 392 parsed_args["rotating"], 393 parsed_args["ocean_vel_fac"], 394 parsed_args["total_hours"], 395 i) 396 Plots.plot!(t/(60.*60.), ice_flux) 397 398 time_elapsed_while_jammed = 0. 399 for it=(length(t) - 1):-1:1 400 if ice_flux[it] ≈ ice_flux[end] 401 time_elapsed_while_jammed = t[end] - t[it] 402 end 403 end 404 405 if time_elapsed_while_jammed > 60.*60. 406 t_jam[i] = t[end] - time_elapsed_while_jammed 407 info("simulation $i jammed at t = $(t_jam[i]/(60.*60.)) h") 408 end 409 410 end 411 Plots.title!(parsed_args["id"]) 412 Plots.xaxis!("Time [h]") 413 Plots.yaxis!("Cumulative ice throughput [kg]") 414 415 Plots.savefig(parsed_args["id"]) 416 Plots.closeall() 417 418 jam_fraction = zeros(length(t)) 419 for it=1:length(t) 420 for i=1:nruns 421 if t_jam[i] <= t[it] 422 jam_fraction[it] += 1./float(nruns) 423 end 424 end 425 end 426 Plots.plot(t/(60.*60.), jam_fraction) 427 Plots.title!(parsed_args["id"] * ", N = " * string(nruns)) 428 Plots.xaxis!("Time [h]") 429 Plots.yaxis!("Fraction of runs jammed [-]") 430 Plots.savefig(parsed_args["id"] * "-jam_fraction.pdf") 431 end 432 433 main()