tsimulation_init.jl - seaice-experiments - sea ice experiments using Granular.jl
(HTM) git clone git://src.adamsgaard.dk/seaice-experiments
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) 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 generation
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 coefficient [-]
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 -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 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.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()