tcompaction.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
       ---
       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 wall
           63         Granular.run!(sim, single_step=true)
           64     end
           65 
           66     #Granular.writeSimulation(sim,filename = "compaction-N$(N)Pa/comp$(filen).jld2")
           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 labels
           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")