deform_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)