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