tpipeline_test.jl - granular-basin - tectonic deformation experiments with Granular.jl
(HTM) git clone git://src.adamsgaard.dk/granular-basin
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
---
tpipeline_test.jl (8164B)
---
1 import Granular
2 import JLD2
3 import PyPlot
4 import Dates
5
6 include("indent.jl")
7
8 #call this script as "@elapsed include("initialization.jl")" in order to time it
9
10 t_start = Dates.now()
11
12 t_init = 0.5 # duration of initialization [s]
13 t_stack = 0.4 #duration for each stacking
14
15 g = [0.,-9.8]
16
17 nx = 25 #125
18 ny = 8 #80
19 stacks = 0
20
21 ngrains = nx*ny*(stacks+1)
22
23 suffix = "_test"
24
25 r_min = 0.05
26 r_max = r_min*sqrt(2) #max grain size is double area of smallest grain size
27
28 SimSettings = Dict()
29
30 SimSettings["nx"] = nx
31 SimSettings["ny"] = ny
32 SimSettings["r_min"] = r_min
33 SimSettings["r_max"] = r_max
34
35 #JLD2.save("SimSettings$(ngrains)$(suffix).jld2", SimSettings)
36
37 gsd_type = "powerlaw"
38 gsd_powerlaw_exponent = -1.8
39 gsd_seed = 3
40
41 # mechanical properties of grains
42 youngs_modulus = 2e7 #elastic modulus
43 poissons_ratio = 0.185 #shear stiffness ratio
44 tensile_strength = 0.0 #strength of bonds between grains
45 contact_dynamic_friction = 0.4 #friction between grains
46 rotating = true #can grains rotate or not
47
48 sim = Granular.createSimulation(id="init")
49 sim.id = "init$(ngrains)$(suffix)"
50
51 Granular.regularPacking!(sim, #simulation object
52 [nx, ny], #number of grains along x and y axis
53 r_min, #smallest grain size
54 r_max, #largest grain size
55 tiling="triangular", #how the grains are tiled
56 origo = [0.0,0.0],
57 size_distribution=gsd_type,
58 size_distribution_parameter=gsd_powerlaw_exponent,
59 seed=gsd_seed)
60
61 for grain in sim.grains #go through all grains in sim
62 grain.youngs_modulus = youngs_modulus
63 grain.poissons_ratio = poissons_ratio
64 grain.tensile_strength = tensile_strength
65 grain.contact_dynamic_friction = contact_dynamic_friction
66 grain.rotating = rotating
67 end
68
69 Granular.fitGridToGrains!(sim, sim.ocean)
70
71 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "north south",
72 verbose=false)
73 #Granular.setGridBoundaryConditions!(sim.ocean, "periodic", "east west")
74 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "east west")
75
76
77 for grain in sim.grains
78 Granular.addBodyForce!(grain, grain.mass*g)
79 Granular.disableOceanDrag!(grain)
80 grain.contact_viscosity_normal = 1e4 # N/(m/s)
81 end
82
83 Granular.setTimeStep!(sim)
84
85 Granular.setTotalTime!(sim, t_init)
86
87 Granular.setOutputFileInterval!(sim, .01)
88
89 Granular.run!(sim)
90
91 #Granular.writeSimulation(sim,
92 # filename = "init$(ngrains)$(suffix).jld2",
93 # folder = "init$(ngrains)$(suffix)")
94
95 #Stack it on top of each other
96
97 temp = deepcopy(sim)
98
99 for i = 1:stacks
100
101 global y_top = -Inf
102 for grain in sim.grains
103 if y_top < grain.lin_pos[2] #+ grain.contact_radius
104 global y_top = grain.lin_pos[2] #+ grain.contact_radius
105 end
106 end
107
108
109 for grain in temp.grains
110
111 lin_pos_lifted = [0.0,0.0]
112
113 lin_pos_lifted[1] = grain.lin_pos[1]
114 lin_pos_lifted[2] = grain.lin_pos[2] + y_top + r_max
115
116
117 Granular.addGrainCylindrical!(sim,
118 lin_pos_lifted,
119 grain.contact_radius,
120 grain.thickness,
121 youngs_modulus = grain.youngs_modulus,
122 poissons_ratio = grain.poissons_ratio,
123 tensile_strength = grain.tensile_strength,
124 contact_dynamic_friction = grain.contact_dynamic_friction,
125 verbose = false)
126
127 end
128
129 Granular.fitGridToGrains!(sim,sim.ocean,verbose=false)
130 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "north south",
131 verbose=false)
132 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "east west",
133 verbose=false)
134 Granular.setTotalTime!(sim,t_stack)
135 Granular.setTimeStep!(sim)
136 Granular.setOutputFileInterval!(sim, .01)
137
138 for grain in sim.grains
139 Granular.addBodyForce!(grain, grain.mass*g)
140 Granular.disableOceanDrag!(grain)
141 end
142
143 #Granular.resetTime!(sim)
144 sim.time_iteration = 0
145 sim.time = 0.0
146 sim.file_time_since_output_file = 0.
147
148 Granular.setTotalTime!(sim, t_stack)
149 Granular.setTimeStep!(sim)
150 Granular.setOutputFileInterval!(sim, .01)
151
152 Granular.run!(sim)
153
154 end
155
156 # add a lower boundary consisting of grains bound together
157 # do this by creating a new simulation where the grains are bonded using
158 # findContacts! after creating overlapping grains. Then set their contact
159 # strengths to an infinite amount, but do not allow new contacts
160
161 carpet = Granular.createSimulation(id="init_carpet")
162
163 bot_r = r_min #radius of bottom layer grains
164
165 left_edge = round(sim.ocean.origo[1],digits=2)
166 length = round(sim.ocean.L[1],digits=2)
167 right_edge = left_edge+length
168
169 for i = left_edge+(bot_r/2):bot_r*1.99:left_edge+length
170
171 bot_pos = [i,round(sim.ocean.origo[2]-bot_r,digits=2)]
172
173 Granular.addGrainCylindrical!(carpet,
174 bot_pos,
175 bot_r,
176 0.1,
177 verbose = false,
178 tensile_strength = Inf,
179 shear_strength = Inf,
180 contact_stiffness_normal = Inf,
181 contact_stiffness_tangential = Inf)
182
183 end
184
185 Granular.findContactsAllToAll!(carpet) #find the grain contacts
186
187 #carpet.grains[10].fixed = true
188 #carpet.grains[10].lin_vel[1:2] = [0.0,0.0]
189
190
191 #add the carpet to the main simulation object
192 append!(sim.grains,carpet.grains)
193
194
195 #fit the ocean to the added grains and run to let the basin settle
196
197 Granular.fitGridToGrains!(sim,sim.ocean,verbose=false)
198
199 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "north south",
200 verbose=false)
201
202 sim.time_iteration = 0
203 sim.time = 0.0
204 sim.file_time_since_output_file = 0.
205 Granular.setTotalTime!(sim, 0.2)
206 Granular.setTimeStep!(sim)
207 Granular.setOutputFileInterval!(sim, .01)
208
209 Granular.run!(sim)
210
211 Granular.writeSimulation(sim,
212 filename = "stacked$(ngrains)$(suffix).jld2")
213
214 ########## Add indenter #############
215
216 temp_indent = createSimulation("id=temp_indent")
217
218
219 left_edge = round(sim.ocean.origo[1],digits=2)
220 length = round(sim.ocean.L[1],digits=2)
221
222 width = length/3
223 hw_ratio = 0.2
224 init_vertex_pos = [(length+left_edge)/2,-0.2]
225 grain_radius = 0.05
226
227
228 vertex_x = init_vertex_pos[1]
229 vertex_y = width*hw_ratio*sin((pi/width)*vertex_x)
230
231
232 for i = 0:grain_radius*2:width#manipulate the ocean grid
233
234 x_pos = i
235
236 y_pos = width*hw_ratio*sin(pi/width*x_pos)
237
238 Granular.addGrainCylindrical!(temp_indent,
239 [x_pos+init_vertex_pos[1]-width/2,y_pos+vertex_y+init_vertex_pos[2]],
240 grain_radius,
241 0.1,
242 fixed = true,
243 lin_vel = [0.0,0.5])
244 end
245
246 append!(sim.grains,temp_indent.grains)
247
248 Granular.fitGridToGrains!(sim,
249 sim.ocean,
250 north_padding = 3.0,
251 verbose=false)
252
253
254 sim.time_iteration = 0
255 sim.time = 0.0
256 sim.file_time_since_output_file = 0.
257 Granular.setTotalTime!(sim, 0.5)
258 Granular.setTimeStep!(sim)
259 Granular.setOutputFileInterval!(sim, .01)
260
261 while sim.time < sim.time_total
262 for grain in carpet.grains
263 if grain.lin_vel[2] < 0 #&& grain.lin_pos[2] < r_min #find some lower boundary here?
264 grain.lin_vel[1:2] = [0.,0.]
265 end
266 end
267 Granular.run!(sim,single_step=true)
268 end
269
270 Granular.writeSimulation(sim,
271 filename = "indented$(ngrains)_test.jld2")
272
273
274 #print time elapsed
275 t_now = Dates.now()
276 dur = Dates.canonicalize(t_now-t_start)
277 print("Time elapsed: ",dur)