compact_basin.jl (3279B)
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 9 id = "simulation40000" # id of simulation to load 10 N = 20e3 # amount of stress to be applied 11 t_comp = 5.0 # compaction max duration [s] 12 13 sim = Granular.readSimulation("$(id)/init.jld2") 14 SimSettings = JLD2.load("$(id)/SimSettings.jld2") 15 #sim = Granular.readSimulation("$(id)/comp.jld2") #use this if continuing already finished compaction 16 17 cd("$id") 18 sim.id = "compaction-N$(N)Pa" 19 SimSettings["N"] = N 20 21 22 Granular.zeroKinematics!(sim) 23 24 y_top = -Inf 25 for grain in sim.grains 26 grain.contact_viscosity_normal = 0 27 if y_top < grain.lin_pos[2] + grain.contact_radius 28 global y_top = grain.lin_pos[2] + grain.contact_radius 29 end 30 end 31 32 #sim.walls = Granular.WallLinearFrictionless[] # remove existing walls 33 34 Granular.addWallLinearFrictionless!(sim, [0., 1.],y_top, 35 bc="normal stress", 36 normal_stress=-N, 37 contact_viscosity_normal=1e3) 38 39 Granular.fitGridToGrains!(sim,sim.ocean) 40 41 y_bot = Inf 42 for grain in sim.grains 43 if y_bot > grain.lin_pos[2] - grain.contact_radius 44 global y_bot = grain.lin_pos[2] - grain.contact_radius 45 end 46 end 47 48 #fixed_thickness = 2. * SimSettings["r_max"] 49 #for grain in sim.grains 50 # if grain.lin_pos[2] <= fixed_thickness 51 # grain.fixed = true # set x and y acceleration to zero 52 # end 53 #end 54 55 Granular.resetTime!(sim) 56 #sim.time_iteration = 0 57 #sim.time = 0.0 58 #sim.file_time_since_output_file = 0. 59 Granular.setTotalTime!(sim,t_comp) 60 61 time = Float64[] 62 compaction = Float64[] 63 effective_normal_stress = Float64[] 64 65 saved5sec = false 66 saved10sec = false 67 68 while sim.time < sim.time_total 69 70 for i = 1:100 # run for a while before measuring the state of the top wall 71 Granular.run!(sim, single_step=true) 72 end 73 74 if sim.time > 5.0 && saved5sec == false 75 Granular.writeSimulation(sim,filename = "comp5sec.jld2") 76 global saved5sec = true 77 end 78 79 if sim.time > 10.0 && saved10sec == false 80 Granular.writeSimulation(sim,filename = "comp10sec.jld2") 81 global saved10sec = true 82 end 83 84 append!(time, sim.time) 85 append!(compaction, sim.walls[1].pos) 86 append!(effective_normal_stress, Granular.getWallNormalStress(sim)) 87 88 end 89 90 defined_normal_stress = ones(size(effective_normal_stress,1)) * 91 Granular.getWallNormalStress(sim, stress_type="effective") 92 93 PyPlot.subplot(211) 94 PyPlot.subplots_adjust(hspace=0.0) 95 ax1 = PyPlot.gca() 96 PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x tick labels 97 PyPlot.plot(time, compaction) 98 PyPlot.ylabel("Top wall height [m]") 99 PyPlot.subplot(212, sharex=ax1) 100 PyPlot.plot(time, defined_normal_stress) 101 PyPlot.plot(time, effective_normal_stress) 102 PyPlot.xlabel("Time [s]") 103 PyPlot.ylabel("Normal stress [Pa]") 104 PyPlot.savefig(sim.id * "-time_vs_compaction-stress.pdf") 105 PyPlot.clf() 106 107 cd("..") 108 109 JLD2.save("simulation$(ngrains)/SimSettings.jld2", SimSettings) 110 111 Granular.writeSimulation(sim,filename = "$(id)/comp.jld2") 112 113 #Granular.writeSimulation(sim,filename = "$(id)/comp10sec.jld2") 114 115 # print time elapsed 116 t_now = Dates.now() 117 dur = Dates.canonicalize(t_now-t_start) 118 print("Time elapsed: ",dur)