tdeform_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
---
tdeform_basin.jl (6110B)
---
1 import Granular
2 import JLD2
3 import Dates
4
5 t_start = Dates.now()
6
7 # User defined settings
8
9 id = "simulation40000" # folder name of simulation
10
11 hw_ratio = 0.2 # height/width ratio of indenter
12 def_time = 4.0 # time spent deforming
13
14 shortening = true # true if walls should be introduced. false if only a diapir
15
16 shortening_type = "fixed" # type of shortening should be "iterative" or "fixed"
17
18 shortening_ratio = 0.05 # ratio of shortening of of basin, if shortening_type
19 # is "fixed". 0.10 would mean basin is shortened by 10%
20
21
22 save_type = "iterative" # "iterative" or "overwrite"
23
24 boomerang_vel = 0.1 # upward velocity of the indeter
25
26
27 t_start = Dates.now()
28
29 sim = Granular.readSimulation("$(id)/layered.jld2")
30 SimSettings = SimSettings = JLD2.load("$(id)/SimSettings.jld2")
31
32 for grain in sim.grains
33 grain.enabled = true
34 grain.fixed = false
35 end
36
37 y_bot_pre = Inf
38 for grain in sim.grains
39 if y_bot_pre > grain.lin_pos[2] - grain.contact_radius
40 global y_bot_pre = grain.lin_pos[2] - grain.contact_radius
41 end
42 end
43
44 # Add Indenter
45 temp_indent = Granular.createSimulation("id=temp_indent")
46
47 left_edge = round(sim.ocean.origo[1],digits=2)
48 length = round(sim.ocean.L[1],digits=2)
49
50 width = length/3
51 init_vertex_pos = [(length+left_edge)/2,-0.2]
52 grain_radius = SimSettings["r_min"]
53
54 vertex_x = init_vertex_pos[1]
55 vertex_y = width*hw_ratio*sin((pi/width)*vertex_x)
56
57
58
59 for i = 0:grain_radius*2:width
60
61 x_pos = i
62
63 y_pos = width*hw_ratio*sin(pi/width*x_pos)
64
65 Granular.addGrainCylindrical!(temp_indent,
66 [x_pos+init_vertex_pos[1]-width/2,y_pos+vertex_y+init_vertex_pos[2]],
67 grain_radius,
68 0.1,
69 fixed = true,
70 lin_vel = [0.0,boomerang_vel],
71 color = -1)
72 end
73
74 append!(sim.grains,temp_indent.grains)
75
76 Granular.fitGridToGrains!(sim,
77 sim.ocean,
78 north_padding = 5.0,
79 verbose=false)
80
81 sim.time_iteration = 0
82 sim.time = 0.0
83 sim.file_time_since_output_file = 0.
84
85 y_bot = Inf
86 for grain in sim.grains
87 if y_bot > grain.lin_pos[2] - grain.contact_radius
88 global y_bot = grain.lin_pos[2] - grain.contact_radius
89 end
90 end
91 Granular.setTotalTime!(sim,def_time)
92 Granular.setTimeStep!(sim)
93 Granular.setOutputFileInterval!(sim, .01)
94 Granular.resetTime!(sim)
95
96 cd("$id")
97 if save_type == "overwrite"
98 sim.id = "deformed"
99 end
100
101 if save_type == "iterative"
102 global save_index = 1
103 while isfile("deformed$(save_index).jld2") == true
104 global save_index += 1
105 end
106 sim.id = "deformed$(save_index)"
107 end
108
109
110 sim.walls = Granular.WallLinearFrictionless[] # remove existing walls
111
112 #find the edge grains of the carpet
113 left_edge = -Inf
114 right_edge = Inf
115 for i = 1:size(sim.grains,1)
116 if left_edge < sim.grains[i].lin_pos[1] + sim.grains[i].contact_radius
117 global left_edge = sim.grains[i].lin_pos[1] + sim.grains[i].contact_radius
118 global left_edge_index = deepcopy(i)
119 end
120 if right_edge > sim.grains[i].lin_pos[1] - sim.grains[i].contact_radius
121 global right_edge = sim.grains[i].lin_pos[1] - sim.grains[i].contact_radius
122 global right_edge_index = deepcopy(i)
123 end
124 end
125
126
127 #add walls to the east and west
128 Granular.addWallLinearFrictionless!(sim,[1.,0.],
129 left_edge,
130 bc = "velocity")
131
132 Granular.addWallLinearFrictionless!(sim,[1.,0.],
133 right_edge,
134 bc = "velocity")
135
136 #add wall beneath the carpet
137
138 Granular.addWallLinearFrictionless!(sim, [0.,1.],
139 y_bot_pre,
140 bc = "fixed")
141
142
143
144
145 global checked_done = false
146
147 #sim.walls[1].vel = -boomerang_vel
148 #sim.walls[2].vel = boomerang_vel
149
150 while sim.time < sim.time_total
151
152 if shortening == true
153 if shortening_type == "iterative"
154 if sim.grains[right_edge_index].lin_vel[1] > boomerang_vel/2 && checked_done == false
155 sim.walls[1].vel = -boomerang_vel
156 sim.walls[2].vel = boomerang_vel
157 global checked_done = true
158 end
159
160 #modulate the speed of the compression walls by the speed of the outer grains
161 if abs(sim.walls[2].vel) > boomerang_vel/2 || abs(sim.walls[2].vel) > abs(sim.grains[right_edge_index].lin_vel[1])*0.98
162 #if boomerang_vel < abs(sim.grains[right_edge_index].lin_vel[1]) || abs(sim.walls[2].vel) > boomerang_vel
163 sim.walls[2].vel *= 0.98
164 end
165 #if abs(sim.walls[2].vel) < abs(sim.grains[right_edge_index].lin_vel[1])*0.90
166 if abs(sim.walls[2].vel) < abs(sim.grains[right_edge_index].lin_vel[1])*0.96
167 sim.walls[2].vel *=1.02
168 end
169
170 #if boomerang_vel < abs(sim.grains[left_edge_index].lin_vel[1]) || abs(sim.walls[1].vel) > boomerang_vel
171 if abs(sim.walls[1].vel) > boomerang_vel/2 || abs(sim.walls[1].vel) > abs(sim.grains[left_edge_index].lin_vel[1])*0.98
172 sim.walls[1].vel *= 0.98
173 end
174 #if abs(sim.walls[1].vel) < abs(sim.grains[left_edge_index].lin_vel[1])*0.90
175 if abs(sim.walls[1].vel) < abs(sim.grains[left_edge_index].lin_vel[1])*0.96
176 sim.walls[1].vel *=1.02
177 end
178 end
179
180 if shortening_type == "fixed" && checked_done == false
181 wall_vel = (length*shortening_ratio)/sim.time_total
182 sim.walls[1].vel = -wall_vel/2
183 sim.walls[2].vel = wall_vel/2
184 global checked_done = true
185 end
186 end
187
188 Granular.run!(sim,single_step = true)
189
190 end
191
192 cd("..")
193
194 Granular.writeSimulation(sim,
195 filename = "$(id)/deformed$(save_index).jld2")
196
197 #print time elapsed
198 t_now = Dates.now()
199 dur = Dates.canonicalize(t_now-t_start)
200 print("Time elapsed: ",dur)