simulation_init.jl (7034B)
1 #/usr/bin/env julia 2 ENV["MPLBACKEND"] = "Agg" 3 import Granular 4 import PyPlot 5 import ArgParse 6 import Random 7 8 const render = false 9 10 function parse_command_line() 11 s = ArgParse.ArgParseSettings() 12 ArgParse.@add_arg_table s begin 13 "--E" 14 help = "Young's modulus [Pa]" 15 arg_type = Float64 16 default = 2e7 17 "--nu" 18 help = "Poisson's ratio [-]" 19 arg_type = Float64 20 default = 0.285 21 "--mu_d" 22 help = "dynamic friction coefficient [-]" 23 arg_type = Float64 24 default = 0.3 25 "--tensile_strength" 26 help = "the maximum tensile strength [Pa]" 27 arg_type = Float64 28 default = 0. 29 "--r_min" 30 help = "minimum grain radius [m]" 31 arg_type = Float64 32 default = 5. 33 "--r_max" 34 help = "maximum grain radius [m]" 35 arg_type = Float64 36 default = 1.35e3 37 default = 50. 38 "--thickness" 39 help = "grain thickness [m]" 40 arg_type = Float64 41 default = 1. 42 "--rotating" 43 help = "allow the grains to rotate" 44 arg_type = Bool 45 default = true 46 "--shearvel" 47 help = "shear velocity [m/s]" 48 arg_type = Float64 49 default = 1.0 50 "--seed" 51 help = "seed for the pseudo RNG" 52 arg_type = Int 53 default = 1 54 "id" 55 help = "simulation id" 56 required = true 57 end 58 return ArgParse.parse_args(s) 59 end 60 61 function report_args(parsed_args) 62 println("Parsed args:") 63 for (arg,val) in parsed_args 64 println(" $arg => $val") 65 end 66 end 67 68 function run_simulation(id::String, 69 E::Float64, 70 nu::Float64, 71 mu_d::Float64, 72 tensile_strength::Float64, 73 r_min::Float64, 74 r_max::Float64, 75 thickness::Float64, 76 rotating::Bool, 77 shearvel::Float64, 78 seed::Int) 79 80 ############################################################################ 81 #### SIMULATION PARAMETERS # 82 ############################################################################ 83 84 # Common simulation identifier 85 id_prefix = id 86 87 # Grain package geometry during initialization 88 nx = 10 # Grains along x (horizontal) 89 ny = 20 # Grains along y (vertical) 90 91 # Grain-size parameters 92 gsd_type = "powerlaw" # "powerlaw" or "uniform" 93 gsd_powerlaw_exponent = -1.8 # GSD power-law exponent 94 gsd_seed = seed # Value to seed random-size generation 95 96 # Grain mechanical properties 97 youngs_modulus = E # Elastic modulus [Pa] 98 poissons_ratio = nu # Shear-stiffness ratio [-] 99 contact_dynamic_friction = mu_d # Coulomb-frictional coefficient [-] 100 101 # Simulation duration of individual steps [s] 102 t_init = 800.0 103 104 # Temporal interval between output files 105 file_dt = 2.0 106 107 ############################################################################ 108 #### Step 1: Create a loose granular assemblage and let it settle at -y # 109 ############################################################################ 110 sim = Granular.createSimulation("$(id_prefix)-init") 111 Granular.removeSimulationFiles(sim) 112 113 #Granular.regularPacking!(sim, [nx, ny], r_min, r_max, 114 #size_distribution=gsd_type, 115 #size_distribution_parameter=gsd_powerlaw_exponent, 116 #seed=gsd_seed) 117 118 # Create a grid for contact searching spanning the extent of the grains 119 sim.ocean = Granular.createRegularOceanGrid([nx, ny, 1], 120 [nx*2*r_max, 121 ny*2*r_max, 122 2*r_max]) 123 124 # Add grains into the initial assemblage 125 Random.seed!(seed) 126 for iy=1:size(sim.ocean.xh, 2) 127 for ix=1:size(sim.ocean.xh, 1) 128 129 for it=1:25 # number of grains to try adding per cell 130 131 r_rand = Granular.randpower(1, gsd_powerlaw_exponent, 132 r_min, r_max) 133 Granular.sortGrainsInGrid!(sim, sim.ocean) 134 pos = Granular.findEmptyPositionInGridCell(sim, sim.ocean, ix, 135 iy, r_rand, 136 n_iter=100, 137 seed=seed*it) 138 if !(typeof(pos) == Array{Float64, 1}) 139 continue 140 end 141 142 Granular.addGrainCylindrical!(sim, pos, r_rand, 1.0, 143 verbose=false) 144 end 145 end 146 end 147 148 # Make the top and bottom boundaries impermeable, and the side boundaries 149 # periodic, which will come in handy during shear 150 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable") 151 Granular.setGridBoundaryConditions!(sim.ocean, "periodic", "east west") 152 153 # Set grain mechanical properties 154 for grain in sim.grains 155 grain.youngs_modulus = youngs_modulus 156 grain.poissons_ratio = poissons_ratio 157 grain.tensile_strength = tensile_strength 158 grain.contact_static_friction = contact_dynamic_friction # mu_s = mu_d 159 grain.contact_dynamic_friction = contact_dynamic_friction 160 grain.rotating = rotating 161 end 162 163 # Drag grains towards -y with the fluid grid 164 sim.ocean.u[:, :, 1, 1] = (rand(nx + 1, ny + 1) .- 0.5).*r_max*0.01 165 sim.ocean.v[:, :, 1, 1] .= -(sim.ocean.yq[end,end]-sim.ocean.xq[1,1])/t_init 166 167 Granular.setTimeStep!(sim) 168 Granular.setTotalTime!(sim, t_init) 169 Granular.setOutputFileInterval!(sim, file_dt) 170 if render 171 Granular.plotGrains(sim, verbose=true) 172 end 173 174 Granular.run!(sim) 175 if render 176 Granular.render(sim) 177 end 178 179 # Create GSD plot to signify that simulation is complete 180 Granular.writeSimulation(sim) 181 if render 182 Granular.plotGrainSizeDistribution(sim) 183 end 184 end 185 186 function main() 187 parsed_args = parse_command_line() 188 report_args(parsed_args) 189 190 seed = parsed_args["seed"] 191 run_simulation(parsed_args["id"] * "-seed" * string(seed), 192 parsed_args["E"], 193 parsed_args["nu"], 194 parsed_args["mu_d"], 195 parsed_args["tensile_strength"], 196 parsed_args["r_min"], 197 parsed_args["r_max"], 198 parsed_args["thickness"], 199 parsed_args["rotating"], 200 parsed_args["shearvel"], 201 seed) 202 end 203 204 main()