tinit-assemblage.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
---
tinit-assemblage.jl (8835B)
---
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 # Common simulation identifier
102 id_prefix = id
103
104 # Grain mechanical properties
105 youngs_modulus = E # Elastic modulus [Pa]
106 poissons_ratio = nu # Shear-stiffness ratio [-]
107 contact_dynamic_friction = mu_d # Coulomb-frictional coefficient [-]
108
109 # Simulation duration [s]
110 t_total = 60.0 * 60.0 * 4.0
111
112 # Temporal interval between output files
113 file_dt = 18.0
114
115 # Misc parameters
116 water_density = 1000.
117
118 # World size [m]
119 L = [5e3, 2e4, 10.0]
120
121 Random.seed!(seed)
122
123 ############################################################################
124 #### Create ice floes
125 ############################################################################
126 sim = Granular.createSimulation("$(id_prefix)")
127 Granular.removeSimulationFiles(sim)
128
129 # Create box edges of ice floes with r_min radius, moving inward
130 r_walls = r_min
131 #= for x in LinRange(0.0, L[1] - r_walls, Int(ceil(L[1]/(r_walls*2)))) =#
132 #= Granular.addGrainCylindrical!(sim, [x, r_min], =#
133 #= r_walls, thickness, fixed=true, =#
134 #= lin_vel=[0.0, compressive_velocity/2.], =#
135 #= verbose=false) =#
136 #= Granular.addGrainCylindrical!(sim, [x, L[2] - r_min], =#
137 #= r_walls, thickness, fixed=true, =#
138 #= lin_vel=[0.0, -compressive_velocity/2.], =#
139 #= verbose=false) =#
140 #= end =#
141 #for y in LinRange(r_walls, L[2] - r_walls, Int(ceil(L[2]/(r_walls*2))))
142 # Granular.addGrainCylindrical!(sim, [r_min, y],
143 # r_walls, thickness, fixed=true,
144 # #lin_vel=[compressive_velocity/2., 0.0],
145 # verbose=false)
146 # Granular.addGrainCylindrical!(sim, [L[1] - r_min, y],
147 # r_walls, thickness, fixed=true,
148 # lin_vel=[-compressive_velocity/2., 0.0],
149 # verbose=false)
150 #end
151
152 # Create a grid for contact searching spanning the extent of the grains
153 n = [Int(ceil(L[1]/r_max/2)), Int(ceil(L[2]/r_max/2)), 2]
154 sim.ocean = Granular.createRegularOceanGrid(n, L)
155 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "-x +x")
156 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "-y +y")
157
158 # Add unconstrained ice floes
159 #= Granular.regularPacking!(sim, =#
160 #= n, =#
161 #= r_min, r_max, =#
162 #= seed=seed) =#
163 Granular.irregularPacking!(sim,
164 radius_max=r_max, radius_min=r_min,
165 thickness=thickness,
166 #binary_radius_search=true,
167 seed=seed)
168
169 # Set grain mechanical properties
170 for grain in sim.grains
171 grain.youngs_modulus = youngs_modulus
172 grain.poissons_ratio = poissons_ratio
173 grain.tensile_strength = abs(tensile_strength +
174 randn()*tensile_strength_std_dev)
175 grain.shear_strength = abs(shear_strength +
176 randn()*tensile_strength_std_dev)
177 grain.contact_static_friction = contact_dynamic_friction # mu_s = mu_d
178 grain.contact_dynamic_friction = contact_dynamic_friction
179 grain.fracture_toughness = fracture_toughness
180 grain.rotating = rotating
181 grain.ocean_drag_coeff_horiz = 0.0 # don't include drag on ends
182 end
183
184 # Create GSD plot to signify that simulation is complete
185 Granular.writeSimulation(sim)
186 render && Granular.plotGrainSizeDistribution(sim)
187
188 # Add a dynamic wall to the top which adds a normal stress downwards.
189 # The normal of this wall is downwards, and we place it at the top of
190 # the granular assemblage. Here, the fluid drag is also disabled
191 y_top = -Inf
192 for grain in sim.grains
193 if y_top < grain.lin_pos[2] + grain.contact_radius
194 y_top = grain.lin_pos[2] + grain.contact_radius
195 end
196 end
197 Granular.addWallLinearFrictionless!(sim, [0., 1.], y_top,
198 bc="normal stress",
199 normal_stress=-20e3)
200 sim.walls[1].mass *= 0.1 # set wall mass to 10% of total grain mass
201 @info("Placing top wall at y=$y_top")
202
203 # Resize the grid to span the current state
204 Granular.fitGridToGrains!(sim, sim.ocean)
205
206 Granular.setTimeStep!(sim)
207 Granular.setTotalTime!(sim, t_total)
208 Granular.setOutputFileInterval!(sim, file_dt)
209 render && Granular.plotGrains(sim, verbose=true)
210
211 ############################################################################
212 #### RUN THE SIMULATION
213 ############################################################################
214
215 Granular.run!(sim)
216 Granular.writeSimulation(sim)
217
218 render && Granular.render(sim, trim=false, animation=true)
219 end
220
221 function main()
222 parsed_args = parse_command_line()
223 report_args(parsed_args)
224
225 seed = parsed_args["seed"]
226 run_simulation(parsed_args["id"] * "-seed" * string(seed),
227 parsed_args["E"],
228 parsed_args["nu"],
229 parsed_args["mu_d"],
230 parsed_args["tensile_strength"],
231 parsed_args["tensile_strength_std_dev"],
232 parsed_args["shear_strength"],
233 parsed_args["fracture_toughness"],
234 parsed_args["r_min"],
235 parsed_args["r_max"],
236 parsed_args["thickness"],
237 parsed_args["rotating"],
238 parsed_args["compressive_velocity"],
239 seed)
240 end
241
242 main()
243 end