 |
tcompaction.jl - granular-basin - tectonic deformation experiments with Granula… |
 |
git clone git://src.adamsgaard.dk/granular-basin |
 |
Log |
|
Files |
|
Refs |
|
README |
|
--- |
|
tcompaction.jl (2506B) |
|
--- |
|
1 import Granular |
|
2 import JLD2 |
|
3 import PyPlot |
|
4 import Dates |
|
5 |
|
6 sim = Granular.readSimulation("stacked60k.jld2") |
|
7 SimSettings = JLD2.load("SimSettings.jld2") |
|
8 |
|
9 N = 20e3 |
|
10 SimSettings["N"] = N |
|
11 |
|
12 JLD2.save("SimSettings.jld2", SimSettings) |
|
13 |
|
14 t_comp = 10.0 #compaction max duration [s] |
|
15 |
|
16 sim.id = "compaction-N$(N)Pa" |
|
17 |
|
18 Granular.zeroKinematics!(sim) |
|
19 |
|
20 y_top = -Inf |
|
21 for grain in sim.grains |
|
22 grain.contact_viscosity_normal = 0 |
|
23 if y_top < grain.lin_pos[2] + grain.contact_radius |
|
24 global y_top = grain.lin_pos[2] + grain.contact_radius |
|
25 end |
|
26 end |
|
27 |
|
28 Granular.addWallLinearFrictionless!(sim, [0., 1.],y_top, |
|
29 bc="normal stress", |
|
30 normal_stress=-N, |
|
31 contact_viscosity_normal=1e3) |
|
32 |
|
33 Granular.fitGridToGrains!(sim,sim.ocean) |
|
34 |
|
35 y_bot = Inf |
|
36 for grain in sim.grains |
|
37 if y_bot > grain.lin_pos[2] - grain.contact_radius |
|
38 global y_bot = grain.lin_pos[2] - grain.contact_radius |
|
39 end |
|
40 end |
|
41 fixed_thickness = 2. * SimSettings["r_max"] |
|
42 for grain in sim.grains |
|
43 if grain.lin_pos[2] <= fixed_thickness |
|
44 grain.fixed = true # set x and y acceleration to zero |
|
45 end |
|
46 end |
|
47 |
|
48 Granular.resetTime!(sim) |
|
49 |
|
50 Granular.setTotalTime!(sim,t_comp) |
|
51 |
|
52 time = Float64[] |
|
53 compaction = Float64[] |
|
54 effective_normal_stress = Float64[] |
|
55 |
|
56 filen = 1 |
|
57 t_start = Dates.now() |
|
58 |
|
59 |
|
60 while sim.time < sim.time_total |
|
61 |
|
62 for i = 1:100 #run for a while before measuring the state of the top… |
|
63 Granular.run!(sim, single_step=true) |
|
64 end |
|
65 |
|
66 #Granular.writeSimulation(sim,filename = "compaction-N$(N)Pa/comp$(f… |
|
67 filen = filen+1 |
|
68 |
|
69 append!(time, sim.time) |
|
70 append!(compaction, sim.walls[1].pos) |
|
71 append!(effective_normal_stress, Granular.getWallNormalStress(sim)) |
|
72 |
|
73 t_now = Dates.now() |
|
74 dur = Dates.canonicalize(t_now-t_start) |
|
75 print("Time elapsed: ",dur) |
|
76 end |
|
77 |
|
78 defined_normal_stress = ones(length(effective_normal_stress)) * |
|
79 Granular.getWallNormalStress(sim, stress_type="effective") |
|
80 |
|
81 PyPlot.subplot(211) |
|
82 PyPlot.subplots_adjust(hspace=0.0) |
|
83 ax1 = PyPlot.gca() |
|
84 PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x tick labe… |
|
85 PyPlot.plot(time, compaction) |
|
86 PyPlot.ylabel("Top wall height [m]") |
|
87 PyPlot.subplot(212, sharex=ax1) |
|
88 PyPlot.plot(time, defined_normal_stress) |
|
89 PyPlot.plot(time, effective_normal_stress) |
|
90 PyPlot.xlabel("Time [s]") |
|
91 PyPlot.ylabel("Normal stress [Pa]") |
|
92 PyPlot.savefig(sim.id * "-time_vs_compaction-stress.pdf") |
|
93 PyPlot.clf() |
|
94 |
|
95 Granular.writeSimulation(sim,filename = "comp.jld2") |
