tinit_basin.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
---
tinit_basin.jl (9058B)
---
1 import Granular
2 import JLD2
3 import PyPlot
4 import Dates
5
6 t_start = Dates.now() # Save the start time, print the end time later.
7
8 ############# Initialization Settings #############
9
10 t_init = 1.5 # duration of initialization [s]
11 t_stack = 1.0 # duration for each stack to settle [s]
12
13 g = [0.,-9.8] # vector for direction and magnitude of gravitational acceleration of grains
14
15 ngrains = 250 # total number of grains
16 aspect_ratio = 4 # should be x times as wide as it is tall
17
18 mkpath("simulation$(ngrains)")
19
20
21 stacks = 2 # number of duplicate stacks on top of the initial grains
22
23 ny = sqrt((ngrains)/aspect_ratio) # number of grain rows
24 nx = aspect_ratio*ny*(stacks+1) # number of grain columns
25
26 ny = Int(round(ny/(stacks+1)))
27 nx = Int(round(nx/(stacks+1)))
28
29
30 r_min = 0.05 # minimum radius of grains
31 r_max = r_min*sqrt(2) # max radius of grains, double the area of minimum
32
33 # Grain-size distribution parameters
34 gsd_type = "powerlaw" # type grain-size distribution
35 gsd_powerlaw_exponent = -1.8 # exponent if powerlaw is used for grain-size distribution
36 gsd_seed = 3 # random seed for the distribution of grain-sizes
37
38 # Initial mechanical properties of grains
39 youngs_modulus = 2e7 # elastic modulus
40 poissons_ratio = 0.185 # shear stiffness ratio
41 tensile_strength = 0.0 # strength of bonds between grains
42 contact_dynamic_friction = 0.4 # friction between grains
43 rotating = true # can grains rotate or not
44
45 # Save some settings in a dictionary
46 SimSettings = Dict()
47 SimSettings["nx"] = nx
48 SimSettings["ny"] = ny
49 SimSettings["stacks"] = stacks
50 SimSettings["ngrains"] = ngrains
51 SimSettings["r_min"] = r_min
52 SimSettings["r_max"] = r_max
53
54 # this section has been moved further up
55
56 ############# Initialize simulation and grains #############
57
58 sim = Granular.createSimulation(id="simulation$(ngrains)")
59
60 Granular.regularPacking!(sim, #simulation object
61 [nx, ny], #number of grains along x and y axis
62 r_min, #smallest grain size
63 r_max, #largest grain size
64 tiling="triangular", #how the grains are tiled
65 origo = [0.0,0.0],
66 size_distribution=gsd_type,
67 size_distribution_parameter=gsd_powerlaw_exponent,
68 seed=gsd_seed,
69 color = 1)
70
71
72 # set the indicated mechanical parameters for all grains
73 for grain in sim.grains
74 grain.youngs_modulus = youngs_modulus
75 grain.poissons_ratio = poissons_ratio
76 grain.tensile_strength = tensile_strength
77 grain.contact_dynamic_friction = contact_dynamic_friction
78 grain.rotating = rotating
79 end
80
81 # fit the ocean grid to the grains
82 Granular.fitGridToGrains!(sim, sim.ocean)
83
84 # set the boundary conditions
85 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "north south",
86 verbose=false)
87 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "east west",
88 verbose=false)
89
90 # add gravity to the grains and remove ocean drag
91 for grain in sim.grains
92 Granular.addBodyForce!(grain, grain.mass*g)
93 Granular.disableOceanDrag!(grain)
94 grain.contact_viscosity_normal = 1e4 # N/(m/s)
95 end
96
97 # run the simulation
98 Granular.setTimeStep!(sim) # set appropriate time steps
99 Granular.setTotalTime!(sim, t_init) # set total time
100 Granular.setOutputFileInterval!(sim, .01) # how often vtu files should be outputted
101
102 cd("simulation$(ngrains)")
103
104 sim.id = "init"
105
106 Granular.run!(sim)
107
108
109
110 ############# Stack the initialized grains #############
111
112 temp = deepcopy(sim)
113
114 for i = 1:stacks
115
116 # find the position of the uppermost grain
117 global y_top = -Inf
118
119 for grain in sim.grains
120 if y_top < grain.lin_pos[2]
121 global y_top = grain.lin_pos[2]
122 end
123 end
124
125 # add duplicate grains above the initialized grains
126 for grain in temp.grains
127
128 lin_pos_lifted = [0.0,0.0] # preallocation of position of a 'lifted' grain
129 lin_pos_lifted[1] = grain.lin_pos[1] # x position of duplicate grain
130 lin_pos_lifted[2] = grain.lin_pos[2] + y_top + r_max # y-position of duplicate grain
131
132
133 Granular.addGrainCylindrical!(sim,
134 lin_pos_lifted,
135 grain.contact_radius,
136 grain.thickness,
137 youngs_modulus = grain.youngs_modulus,
138 poissons_ratio = grain.poissons_ratio,
139 tensile_strength = grain.tensile_strength,
140 contact_dynamic_friction = grain.contact_dynamic_friction,
141 verbose = false)
142
143 end
144
145 Granular.fitGridToGrains!(sim,sim.ocean,verbose=false)
146 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "north south",
147 verbose=false)
148 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "east west",
149 verbose=false)
150 Granular.setTotalTime!(sim,t_stack)
151 Granular.setTimeStep!(sim)
152 Granular.setOutputFileInterval!(sim, .01)
153
154 for grain in sim.grains
155 Granular.addBodyForce!(grain, grain.mass*g)
156 Granular.disableOceanDrag!(grain)
157 end
158
159
160 # instead of Granular.resetTime!(sim), the time parameters are reset manually, in order
161 # to allow for the continuation of vtu-file index from earlier
162 sim.time_iteration = 0
163 sim.time = 0.0
164 sim.file_time_since_output_file = 0.
165
166 Granular.setTotalTime!(sim, t_stack)
167 Granular.setTimeStep!(sim)
168 Granular.setOutputFileInterval!(sim, .01)
169
170 Granular.run!(sim) # let the duplicated grains settle.
171
172 end
173
174
175
176
177 ############# Lay a carpet #############
178
179 carpet = Granular.createSimulation(id="init_carpet") # new simulation object for the carpet
180
181 bot_r = r_min # radius of carpet grains
182
183 left_edge = round(sim.ocean.origo[1],digits=2) # west edge of the carpet
184 length = round(sim.ocean.L[1],digits=2) # width of the carpet
185 right_edge = left_edge+length # east edge of the carpet
186
187
188 # Now loop over the carpet grain positions, the loop will create grains that overlap slightly
189 # in order to create the bonds needed
190 # color = 0 is used as a flag for the grains in the carpet
191 for i = left_edge+(bot_r/2):bot_r*1.999:left_edge+length
192
193 bot_pos = [i,round(sim.ocean.origo[2]-bot_r,digits=2)] # position of grain
194
195 Granular.addGrainCylindrical!(carpet,
196 bot_pos,
197 bot_r,
198 0.1,
199 verbose = false,
200 tensile_strength = Inf,
201 shear_strength = Inf,
202 #contact_stiffness_normal = Inf,
203 #contact_stiffness_tangential = Inf,
204 fixed = false,
205 color = 0)
206 end
207
208 #Granular.fitGridToGrains!(carpet,carpet.ocean,verbose=false)
209
210
211
212 Granular.findContactsAllToAll!(carpet) # find the grain contacts
213
214
215
216 append!(sim.grains,carpet.grains) # add the carpet grains to the main simulation object
217 # since the assignment will point to the carpet object, changes made to the carpet
218 # object will appear in the main simulation object
219
220
221 """
222 #reset the grain contacts and make them very old
223
224 for grain in sim.grains
225 grain.contacts[:] .= 0
226 grain.n_contacts = 0
227 end
228
229
230 for grain in sim.grains
231 for ic=1:size(grain.contact_age,1)
232 grain.contact_age[ic] = 1e16
233 end
234 grain.strength_heal_rate = 1 # new bond stengthening
235 end
236 """
237
238 Granular.fitGridToGrains!(sim,sim.ocean,verbose=false) # fit the ocean to the added grains
239
240 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "north south",
241 verbose=false)
242 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "east west",
243 verbose=false)
244
245 #Granular.findContacts!(sim,method="ocean grid")
246
247 # run the simulation shortly, to let the stacked grains settle on the carpet
248 sim.time_iteration = 0
249 sim.time = 0.0
250 sim.file_time_since_output_file = 0.
251 Granular.setTotalTime!(sim, 0.5)
252 Granular.setTimeStep!(sim)
253 Granular.setOutputFileInterval!(sim, .01)
254
255
256 Granular.run!(sim)
257
258
259 # save the simulation and the carpet objects
260
261 cd("..")
262
263 Granular.writeSimulation(sim,
264 filename = "simulation$(ngrains)/init.jld2")
265
266 JLD2.save("simulation$(ngrains)/SimSettings.jld2", SimSettings)
267
268 # print time elapsed
269 t_now = Dates.now()
270 dur = Dates.canonicalize(t_now-t_start)
271 print("Time elapsed: ",dur)