GopherProxy

	tsimulation_mohr-coulomb_plot.jl - seaice-experiments - sea ice experiments usi…
	git clone git://src.adamsgaard.dk/seaice-experiments
	Log
	Files
	Refs
	README
	LICENSE
	---
	tsimulation_mohr-coulomb_plot.jl (4045B)
	---
	1 #/usr/bin/env julia
	2 ENV["MPLBACKEND"] = "Agg"
	3 import PyPlot
	4 import LsqFit
	5
	6 const id_prefix = "mohr_coulomb"
	7 const vel_shear = 1.0
	8
	9 # Normal stresses for consolidation and shear [Pa]
	10 #const N_list = [20e3, 40e3, 80e3, 160e3]
	11 const N_list = [5e3, 10e3, 15e3, 20e3, 30e3, 40e3, 80e3, 160e3]
	12
	13 #time = Float64[]
	14 #shear_strain = Float64[]
	15 shear_stress = Float64[]
	16 #dilation = Float64[]
	17 tau_p = Float64[]
	18 tau_u = Float64[]
	19 sim_id = ""
	20
	21 for kind in ["mu0.3_sigma_c0kPa", "mu0.0_sigma_c200kPa"]
	22 for N in N_list
	23
	24 #mohr_coulomb_mu0.3_sigma_c0kPa.pdf-seed1-shear-N20000.0Pa-vel_s…
	25 sim_id = "$(id_prefix)_$(kind).pdf-seed1-shear-N$(N)Pa-vel_shear…
	26 info(sim_id)
	27 data = readdlm(sim_id * "-data.txt")
	28
	29 #time = data[1,:]
	30 #shear_strain = data[2,:]
	31 shear_stress = data[3,:]
	32 #dilation = data[4,:]
	33
	34 append!(tau_p, maximum(shear_stress[1:100]))
	35 append!(tau_u, mean(shear_stress[end-100:end]))
	36 end
	37 end
	38
	39 # Plot normal stress vs. peak shear friction
	40 PyPlot.figure(figsize=(5, 4))
	41 PyPlot.plot(N_list, tau_p[1:length(N_list)]./N_list,
	42 "xk",
	43 label="Frictional DEM")
	44 PyPlot.plot(N_list, tau_p[length(N_list)+1:end]./N_list,
	45 "+b",
	46 label="Cohesive DEM")
	47 PyPlot.ylabel("Peak shear friction \$\\mu_p\$ [-]")
	48
	49 PyPlot.xlabel("Normal stress \$N\$ [Pa]")
	50 PyPlot.xlim(0., maximum(N_list)*1.1)
	51 PyPlot.legend()
	52 PyPlot.savefig(id_prefix * "-mu_p.pdf")
	53 PyPlot.clf()
	54
	55 # Plot normal stress vs. effective shear friction
	56 PyPlot.plot(N_list./1e3, tau_u[1:length(N_list)]./N_list,
	57 "xk",
	58 label="Frictional DEM")
	59 PyPlot.plot(N_list./1e3, tau_u[length(N_list)+1:end]./N_list,
	60 "+b",
	61 label="Cohesive DEM")
	62 PyPlot.xlabel("Normal stress \$N\$ [kPa]")
	63 PyPlot.ylabel("Effective shear friction \$\\tau_u/N\$ [-]")
	64 #PyPlot.xlim(0., maximum(N_list)*1.1/1e3)
	65 #PyPlot.xlim(minimum(N_list)0.9/1e3, maximum(N_list)1.1/1e3)
	66 PyPlot.legend()
	67 PyPlot.title("(b)")
	68 PyPlot.savefig(id_prefix * "-mu_eff.pdf")
	69 PyPlot.clf()
	70
	71 ### Plot normal stress vs. ultimate shear stress
	72 mohr_coulomb(N, param) = param[2] + param[1]*N
	73
	74 # Frictional DEM
	75 fit = LsqFit.curve_fit(mohr_coulomb, N_list, tau_u[1:length(N_list)], [0…
	76 println("### Frictional DEM ###")
	77 println(" μ_u = $(fit.param[1]), C = $(fit.param[2]/1e3) kPa")
	78 errors = LsqFit.estimate_errors(fit, 0.95)
	79 println(" 95% CI errors:")
	80 println(" μ_u: ± $(errors[1])")
	81 println(" C: ± $(errors[2]/1e3) kPa")
	82 label = @sprintf("\$\\tau_u\$(\$N\$, \$\\mu_u\$ = %.2g±%.2g, \$C\$ = %.…
	83 fit.param[1], errors[1],
	84 fit.param[2]/1e3, errors[2]/1e3)
	85 PyPlot.plot(N_list./1e3, tau_u[1:length(N_list)]./1e3,
	86 "xk",
	87 label="Frictional DEM")
	88 PyPlot.plot(linspace(N_list[1], N_list[end])./1e3,
	89 mohr_coulomb(linspace(N_list[1], N_list[end]), fit.param)./1…
	90 "-k", alpha=0.5, label=label)
	91
	92 # Cohesive DEM
	93 fit = LsqFit.curve_fit(mohr_coulomb, N_list, tau_u[length(N_list)+1:end]…
	94 println("### Cohesive DEM ###")
	95 println(" μ_u = $(fit.param[1]), C = $(fit.param[2]/1e3) kPa")
	96 errors = LsqFit.estimate_errors(fit, 0.95)
	97 println(" 95% CI errors:")
	98 println(" μ_u: ± $(errors[1])")
	99 println(" C: ± $(errors[2]/1e3) kPa")
	100 label = @sprintf("\$\\tau_u\$(\$N\$, \$\\mu_u\$ = %.2g±%.2g, \$C\$ = %.…
	101 fit.param[1], errors[1],
	102 fit.param[2]/1e3, errors[2]/1e3)
	103 PyPlot.plot(N_list./1e3, tau_u[length(N_list)+1:end]./1e3,
	104 "+b",
	105 label="Cohesive DEM")
	106 PyPlot.plot(linspace(N_list[1], N_list[end])./1e3,
	107 mohr_coulomb(linspace(N_list[1], N_list[end]), fit.param)./1…
	108 "--b", alpha=0.5, label=label)
	109
	110 # Annotation
	111 PyPlot.xlabel("Normal stress \$N\$ [kPa]")
	112 PyPlot.ylabel("Shear stress \$\\tau_u\$ [kPa]")
	113 #PyPlot.xlim(minimum(N_list)0.9/1e3, maximum(N_list)1.1/1e3)
	114 PyPlot.legend()
	115 PyPlot.title("(a)")
	116 PyPlot.savefig(id_prefix * "-tau_u.pdf")
	117 PyPlot.clf()