tcompress.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
---
tcompress.jl (8221B)
---
1 #/usr/bin/env julia
2 ENV["MPLBACKEND"] = "Agg"
3 import Granular
4 import PyPlot
5 import ArgParse
6 using Random
7 import DelimitedFiles
8
9 let
10 render = false
11
12 function parse_command_line()
13 s = ArgParse.ArgParseSettings()
14 ArgParse.@add_arg_table s begin
15 "--E"
16 help = "Young's modulus [Pa]"
17 arg_type = Float64
18 default = 2e7
19 "--nu"
20 help = "Poisson's ratio [-]"
21 arg_type = Float64
22 default = 0.285
23 "--mu_d"
24 help = "dynamic friction coefficient [-]"
25 arg_type = Float64
26 default = 0.3
27 "--tensile_strength"
28 help = "the maximum tensile strength [Pa]"
29 arg_type = Float64
30 default = 400e3
31 "--tensile_strength_std_dev"
32 help = "standard deviation of the maximum tensile strength [Pa]"
33 arg_type = Float64
34 default = 0.0
35 "--shear_strength"
36 help = "the maximum shear strength [Pa]"
37 arg_type = Float64
38 default = 200e3
39 "--fracture_toughness"
40 help = "fracture toughness [Pa*m^{1/2}]"
41 arg_type = Float64
42 default = 1e20
43 #default = 1285e3
44 "--r_min"
45 help = "min. grain radius [m]"
46 arg_type = Float64
47 default = 20.0
48 "--r_max"
49 help = "max. grain radius [m]"
50 arg_type = Float64
51 default = 80.0
52 "--thickness"
53 help = "grain thickness [m]"
54 arg_type = Float64
55 default = 1.
56 "--rotating"
57 help = "allow the grains to rotate"
58 arg_type = Bool
59 default = true
60 "--compressive_velocity"
61 help = "compressive velocity [m/s]"
62 arg_type = Float64
63 default = 0.1
64 "--seed"
65 help = "seed for the pseudo RNG"
66 arg_type = Int
67 default = 1
68 "id"
69 help = "simulation id"
70 required = true
71 end
72 return ArgParse.parse_args(s)
73 end
74
75 function report_args(parsed_args)
76 println("Parsed args:")
77 for (arg,val) in parsed_args
78 println(" $arg => $val")
79 end
80 end
81
82 function run_simulation(id::String,
83 E::Float64,
84 nu::Float64,
85 mu_d::Float64,
86 tensile_strength::Float64,
87 tensile_strength_std_dev::Float64,
88 shear_strength::Float64,
89 fracture_toughness::Float64,
90 r_min::Float64,
91 r_max::Float64,
92 thickness::Float64,
93 rotating::Bool,
94 compressive_velocity::Float64,
95 seed::Int)
96
97 ############################################################################
98 #### SIMULATION PARAMETERS
99 ############################################################################
100
101 # Grain mechanical properties
102 youngs_modulus = E # Elastic modulus [Pa]
103 poissons_ratio = nu # Shear-stiffness ratio [-]
104 contact_dynamic_friction = mu_d # Coulomb-frictional coefficient [-]
105
106 # Simulation duration [s]
107 t_total = 60.0 * 60.0 * 2.0
108
109 # Temporal interval between output files
110 file_dt = 18.0
111
112 ############################################################################
113 #### Read prior state
114 ############################################################################
115 sim = Granular.readSimulation("rotating01-seed1/rotating01-seed1.681.jld2")
116 sim.id = id
117 Granular.resetTime!(sim)
118 Granular.zeroKinematics!(sim)
119 Granular.removeSimulationFiles(sim)
120
121 # Set grain mechanical properties
122 for grain in sim.grains
123 grain.youngs_modulus = youngs_modulus
124 grain.poissons_ratio = poissons_ratio
125 grain.tensile_strength = abs(tensile_strength +
126 randn()*tensile_strength_std_dev)
127 grain.shear_strength = abs(shear_strength +
128 randn()*tensile_strength_std_dev)
129 grain.contact_static_friction = contact_dynamic_friction # mu_s = mu_d
130 grain.contact_dynamic_friction = contact_dynamic_friction
131 grain.fracture_toughness = fracture_toughness
132 grain.rotating = rotating
133 grain.ocean_drag_coeff_horiz = 0.0 # don't include drag on ends
134 end
135
136 x_max = -Inf
137 x_min = +Inf
138 y_max = -Inf
139 y_min = +Inf
140 r_max = 0.0
141 for grain in sim.grains
142 r = grain.contact_radius
143
144 if x_min > grain.lin_pos[1] - r
145 x_min = grain.lin_pos[1] - r
146 end
147
148 if x_max < grain.lin_pos[1] + r
149 x_max = grain.lin_pos[1] + r
150 end
151
152 if y_min > grain.lin_pos[2] - r
153 y_min = grain.lin_pos[2] - r
154 end
155
156 if y_max < grain.lin_pos[2] + r
157 y_max = grain.lin_pos[2] + r
158 end
159
160 if r_max < r
161 r_max = r
162 end
163 end
164
165 dx = r_max * 2.0
166 L::Vector{Float64} = [(x_max - x_min)*3.0, y_max - y_min, dx]
167 n = convert(Vector{Int}, floor.(L./dx))
168 sim.ocean = Granular.createRegularOceanGrid(n, L,
169 origo=[x_min - L[1]/3.0, y_min],
170 time=[0.], name="resized")
171 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "-x +x")
172 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "-y +y")
173
174 sim.walls[1].bc = "velocity"
175 sim.walls[1].vel = -(y_max - y_min) * 0.2 / t_total
176
177 # fix grains adjacent to walls to create no-slip BC
178 fixed_thickness = 2.0 * r_max
179 for grain in sim.grains
180 if grain.lin_pos[2] <= y_min + fixed_thickness
181 grain.fixed = true
182 grain.color = 1
183 elseif grain.lin_pos[2] >= y_max - fixed_thickness
184 grain.allow_y_acc = true
185 grain.fixed = true
186 grain.color = 1
187 end
188 end
189
190 Granular.setTimeStep!(sim)
191 Granular.setTotalTime!(sim, t_total)
192 Granular.setOutputFileInterval!(sim, file_dt)
193 render && Granular.plotGrains(sim, verbose=true)
194
195 ############################################################################
196 #### RUN THE SIMULATION
197 ############################################################################
198
199 # Run the consolidation experiment, and monitor top wall position over
200 # time
201 time = Float64[]
202 compaction = Float64[]
203 effective_normal_stress = Float64[]
204 while sim.time < sim.time_total
205
206 for i=1:100
207 Granular.run!(sim, single_step=true)
208 end
209
210 append!(time, sim.time)
211 append!(compaction, sim.walls[1].pos)
212 append!(effective_normal_stress, Granular.getWallNormalStress(sim))
213
214 end
215 Granular.writeSimulation(sim)
216
217 defined_normal_stress = ones(length(effective_normal_stress)) *
218 Granular.getWallNormalStress(sim, stress_type="effective")
219 PyPlot.subplot(211)
220 PyPlot.subplots_adjust(hspace=0.0)
221 ax1 = PyPlot.gca()
222 PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x labels
223 PyPlot.plot(time, compaction)
224 PyPlot.ylabel("Top wall height [m]")
225 PyPlot.subplot(212, sharex=ax1)
226 PyPlot.plot(time, defined_normal_stress)
227 PyPlot.plot(time, effective_normal_stress)
228 PyPlot.xlabel("Time [s]")
229 PyPlot.ylabel("Normal stress [Pa]")
230 PyPlot.savefig(sim.id * "-time_vs_compaction-stress.pdf")
231 PyPlot.clf()
232 DelimitedFiles.writedlm(sim.id * "-data.txt", [time, compaction,
233 effective_normal_stress,
234 defined_normal_stress])
235
236 render && Granular.render(sim, trim=false, animation=true)
237 end
238
239 function main()
240 parsed_args = parse_command_line()
241 report_args(parsed_args)
242
243 seed = parsed_args["seed"]
244 run_simulation(parsed_args["id"] * "-seed" * string(seed),
245 parsed_args["E"],
246 parsed_args["nu"],
247 parsed_args["mu_d"],
248 parsed_args["tensile_strength"],
249 parsed_args["tensile_strength_std_dev"],
250 parsed_args["shear_strength"],
251 parsed_args["fracture_toughness"],
252 parsed_args["r_min"],
253 parsed_args["r_max"],
254 parsed_args["thickness"],
255 parsed_args["rotating"],
256 parsed_args["compressive_velocity"],
257 seed)
258 end
259
260 main()
261 end