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