seaice-experiments

sea ice experiments using Granular.jl
git clone git://src.adamsgaard.dk/seaice-experiments # fast
git clone https://src.adamsgaard.dk/seaice-experiments.git # slow
Log | Files | Refs | README | LICENSE Back to index

ridging_bulk_plots.jl (1722B)


      1 #/usr/bin/env julia
      2 ENV["MPLBACKEND"] = "Agg"
      3 import PyPlot
      4 using DelimitedFiles
      5 
      6 id_prefix = "ridging_bulk_elastic_plastic1"
      7 #compressive_velocity = [0.2, 0.1, 0.05]
      8 compressive_velocity = [0.1]
      9 fracture_toughness = ["1285e2", "2570e2",
     10                       "1285e3", "2570e3",
     11                       "1285e4", "2570e4"]
     12 fracture_toughness_labels = ["1.285e5", "2.570e5",
     13                              "1.285e6", "2.570e6",
     14                              "1.285e7", "2.570e7"]
     15 
     16 function readTimeSeries(id::String,
     17                         cv::Float64)
     18 
     19     data = readdlm(id * "-data.txt")  # file is: time, N, τ
     20 
     21     return data[1,:], data[1,:].*cv/800.0, data[2,:]./1e3, data[3,:]./1e3
     22 end
     23 
     24 PyPlot.figure(figsize=(4,4))
     25 i = 6
     26 for cv in compressive_velocity
     27     for K in fracture_toughness[end:-1:1]
     28         t, γ, N, τ = readTimeSeries(id_prefix * "-K$(K)-cv$(cv)-seed1",
     29                                     cv)
     30 
     31         #= PyPlot.plot(γ, τ./N, =#
     32         #=             linewidth=0.2, alpha=0.5, =#
     33         #=             label="$K = $K$ Pa m$^{1/2}$") =#
     34         #PyPlot.plot(γ, N,
     35         PyPlot.semilogy(γ, abs.(N),
     36                     linewidth=0.5, alpha=0.75,
     37                     label="\$K_\\mathrm{Ic}\$ = $(fracture_toughness_labels[i]) Pa m\$^{1/2}\$")
     38         global i -= 1
     39     end
     40 end
     41 PyPlot.legend()
     42 #= PyPlot.xlabel("Shear strain, \$\\gamma\$ [-]") =#
     43 PyPlot.xlabel("Compressive strain, \$\\epsilon_c\$ [-]")
     44 #= PyPlot.ylabel("Shear friction, $\\mu = \\tau/N$, [-]") =#
     45 PyPlot.ylabel("Bulk compressive stress, \$N\$ [kPa]")
     46 PyPlot.xlim([0.0, 0.45])
     47 PyPlot.ylim([1e-4, 2e2])
     48 PyPlot.tight_layout()
     49 PyPlot.savefig("ridging_bulk_uniaxial_strain_N.pdf")
     50 PyPlot.savefig("ridging_bulk_uniaxial_strain_N.png")