GopherProxy

	tcompact_basin.jl - granular-basin - tectonic deformation experiments with Gran…
	git clone git://src.adamsgaard.dk/granular-basin
	Log
	Files
	Refs
	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 continuin…
	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 to…
	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 labe…
	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)