GopherProxy

	tsimulation_init.jl - seaice-experiments - sea ice experiments using Granular.jl
	git clone git://src.adamsgaard.dk/seaice-experiments
	Log
	Files
	Refs
	README
	LICENSE
	---
	tsimulation_init.jl (7087B)
	---
	1 #/usr/bin/env julia
	2 ENV["MPLBACKEND"] = "Agg"
	3 import Granular
	4 import JLD
	5 import PyPlot
	6 import ArgParse
	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 const id_prefix = id
	86
	87 # Grain package geometry during initialization
	88 const nx = 10 # Grains along x (horizontal)
	89 const ny = 20 # Grains along y (vertical)
	90
	91 # Grain-size parameters
	92 const gsd_type = "powerlaw" # "powerlaw" or "uniform"
	93 const gsd_powerlaw_exponent = -1.8 # GSD power-law exponent
	94 const gsd_seed = seed # Value to seed random-size ge…
	95
	96 # Grain mechanical properties
	97 const youngs_modulus = E # Elastic modulus [Pa]
	98 const poissons_ratio = nu # Shear-stiffness ratio [-]
	99 const contact_dynamic_friction = mu_d # Coulomb-frictional coefficie…
	100
	101 # Simulation duration of individual steps [s]
	102 const t_init = 800.0
	103
	104 # Temporal interval between output files
	105 const file_dt = 2.0
	106
	107 ####################################################################…
	108 #### Step 1: Create a loose granular assemblage and let it settle at…
	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_e…
	116 #seed=gsd_seed)
	117
	118 # Create a grid for contact searching spanning the extent of the gra…
	119 sim.ocean = Granular.createRegularOceanGrid([nx, ny, 1],
	120 [nx2r_max,
	121 ny2r_max,
	122 2*r_max])
	123
	124 # Add grains into the initial assemblage
	125 srand(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.ocea…
	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 bound…
	149 # periodic, which will come in handy during shear
	150 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable")
	151 Granular.setGridBoundaryConditions!(sim.ocean, "periodic", "east wes…
	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…
	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_max0.01
	165 sim.ocean.v[:, :, 1, 1] = -(sim.ocean.yq[end,end]-sim.ocean.xq[1,1])…
	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()