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

tau_plots.jl (6162B)


      1 #!/usr/bin/env julia
      2 import ArgParse
      3 ENV["MPLBACKEND"] = "Agg"
      4 using PyPlot
      5 
      6 function plotTau(xvalues::Vector, tau::Vector, tau_s::Vector,
      7                  x_label::String, sim_id::String;
      8                  x_lim::Vector=[-1., -1.],
      9                  y_lim::Vector=[-1., -1.],
     10                  ycenter::Float64=-1,
     11                  shadex::Vector=[-1., -1.],
     12                  xvalues2::Vector=[-1.],
     13                  tau2::Vector=[-1.],
     14                  tau2_s::Vector=[-1.],
     15                  x_label2::String="",
     16                  x_lim2::Vector=[-1., -1.],
     17                  logx::Bool=false,
     18                  cohesive::Bool=false)
     19 
     20     fig = figure(figsize=(5,4))
     21     ax1 = gca()
     22     if length(xvalues2) > 1
     23         ax2 = ax1[:twiny]()
     24 
     25         #new_position = [0.06;0.06;0.77;0.91] # Position Method 2
     26         #ax2[:set_position](top=0.85)
     27     end
     28 
     29     if length(tau2) > 1
     30         #semilogy(xvalues, tau, "xk", label="Frictional DEM")
     31         #semilogy(xvalues, tau2, "+k", label="Cohesive DEM")
     32         ax1[:errorbar](xvalues, tau, yerr=tau_s, fmt="xk",
     33                         label="Frictional DEM")
     34 
     35         if length(xvalues2) > 1
     36             ax2[:errorbar](xvalues2, tau2, yerr=tau2_s, fmt="+b",
     37                            label="Cohesive DEM")
     38             ax1[:legend](fancybox="true", loc="lower right")
     39             ax2[:legend](fancybox="true", loc="upper right")
     40         else
     41             ax1[:errorbar](xvalues, tau2, yerr=tau2_s, fmt="+b",
     42                            label="Cohesive DEM")
     43             ax1[:legend](fancybox="true")
     44         end
     45     else
     46         #semilogy(xvalues, tau, "ok")
     47         if cohesive
     48             errorbar(xvalues, tau, yerr=tau_s, fmt="+b")
     49         else
     50             errorbar(xvalues, tau, yerr=tau_s, fmt="xk")
     51         end
     52     end
     53 
     54     yscale("log")
     55     if logx
     56         xscale("log")
     57     end
     58 
     59     if shadex[2] > 0.
     60         ax1[:axvspan](shadex[1], shadex[2], color="gray", alpha=0.5, lw=0)
     61 
     62         ax1[:text]((shadex[2] - shadex[1])/2. + shadex[1],
     63                     ycenter,
     64                     "No jamming observed",
     65                     horizontalalignment="center",
     66                     verticalalignment="center",
     67                     rotation="vertical")
     68     end
     69 
     70     for x in xvalues
     71         ax1[:plot]([x, x], y_lim[1].*[0.9, 1.1], "-r", lw=4)
     72     end
     73     if length(xvalues2) > 1
     74         for x in xvalues2
     75             ax2[:plot]([x, x], y_lim[2].*[0.9, 1.1], "-r", lw=4)
     76         end
     77     end
     78 
     79     #title(parsed_args["id"])
     80     ax1[:set_xlabel](x_label)
     81     ax1[:set_ylabel]("Time-scale for jamming \$T\$ [h]")
     82     if x_lim[2] > 0.
     83         ax1[:set_xlim](x_lim)
     84     end
     85     if y_lim[2] > 0.
     86         ax1[:set_ylim](y_lim)
     87     end
     88 
     89     if length(xvalues2) > 1
     90         font1 = Dict("color"=>"blue")
     91         ax2[:set_xlabel](x_label2, fontdict=font1)
     92         setp(ax2[:get_xticklabels](), color="blue")
     93 
     94         if x_lim2[2] > 0.
     95             ax2[:set_xlim](x_lim2)
     96         end
     97     end
     98 
     99     tight_layout()
    100 
    101     if length(xvalues2) > 1
    102         savefig(sim_id * "_combined_tau.png")
    103         savefig(sim_id * "_combined_tau.pdf")
    104     else
    105         savefig(sim_id * "_$(x_label)_" * "_tau.png")
    106         savefig(sim_id * "_$(x_label)_" * "_tau.pdf")
    107     end
    108 end
    109 
    110 # mu_sigmac
    111 
    112 data = readdlm("mu_sigmac/mu_sigmac_fits.txt")
    113 tau = data[:,2]
    114 s = data[:,3]
    115 plotTau([0.00, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45],
    116     tau[1:9],
    117     s[1:9],
    118     "Coulomb-friction coefficient \$\\mu\$ [-]", "mu_sigmac",
    119     x_lim=[0.0, 0.5],
    120     y_lim=[1e0, 2e2],
    121     ycenter=1.05e1,
    122     shadex=[0.0, 0.08])
    123 
    124 plotTau(
    125     [    0.,    10.,  100.,  200., 300., 400., 500., 600., 700., 800.],
    126     vcat(tau[1], tau[10:end]),
    127     vcat(s[1], s[10:end]),
    128     "Tensile strength \$\\sigma_c\$ [kPa]", "mu_sigmac",
    129     x_lim=[0, 1000.],
    130     y_lim=[1e-1, 1e2],
    131     ycenter=3e0,
    132     shadex=[0., 150.],
    133     cohesive=true)
    134 
    135 
    136 ## combined mu_sigmac plot
    137 plotTau([0.00, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45],
    138     tau[1:9],
    139     s[1:9],
    140     "Coulomb-friction coefficient \$\\mu\$ [-]",
    141     "mu_sigmac",
    142     #x_lim=[0.0, 0.5],
    143     x_lim=[0.0, 0.9],
    144     y_lim=[1e-1, 2e2],
    145     ycenter=1.05e1,
    146     shadex=[0.0, 0.09],
    147     xvalues2=[    0.,    10.,  100.,  200., 300., 400., 500., 600., 700., 800.],
    148     x_lim2=[50.0, 750.],
    149     tau2=vcat(tau[1], tau[10:end]),
    150     tau2_s=vcat(s[1], s[10:end]),
    151     x_label2="Tensile strength \$\\sigma_c\$ [kPa]")
    152 
    153 
    154 # psd_width
    155 data = readdlm("psd_width/psd_width_fits.txt")
    156 tau = data[:,2]
    157 s = data[:,3]
    158 plotTau(
    159     [0., 1.15e3-8e2, 1.35e3-6e2, 1.55e3-4e2, 1.75e3-2e2]*2.,
    160     tau[5:-1:1],
    161     s[5:-1:1],
    162     "PSD width \$d_{max} - d_{min}\$ [m]", "psd_width",
    163     x_lim=[0, 3300.],
    164     y_lim=[1e-0, 2e2],
    165     ycenter=1e1,
    166     shadex=[0., 500.],
    167     tau2=tau[end:-1:6],
    168     tau2_s=s[end:-1:6])
    169 
    170 # stiffness
    171 #=data = readdlm("stiffness/stiffness_fits.txt")
    172 tau = data[:,2]
    173 s = data[:,3]
    174 plotTau([2e5, 2e6, 2e7, 2e8, 2e9],
    175     tau[1:5],
    176     s[1:5],
    177     "Young's modulus \$E\$ [Pa]", "stiffness",
    178     x_lim=[1e5, 1e10],
    179     y_lim=[1e-1, 2e3],
    180     ycenter=3e0,
    181     shadex=[0., 0.],
    182     tau2=tau[6:end],
    183     tau2_s=s[6:end],
    184     logx=true)=#
    185 
    186 # thickness
    187 data = readdlm("thickness/thickness_fits.txt")
    188 tau = data[:,2]
    189 s = data[:,3]
    190 plotTau([1.0, 1.29, 1.67, 2.15, 2.78, 3.59, 4.64, 5.99, 7.74, 10.],
    191     tau[1:10],
    192     s[1:10],
    193     "Ice-floe thickness \$h\$ [m]", "thickness",
    194     x_lim=[0.9, 11.],
    195     y_lim=[1e0, 1e3],
    196     ycenter=3e0,
    197     shadex=[0.1, 0.11],
    198     tau2=tau[11:end],
    199     tau2_s=s[11:end],
    200     logx=true)
    201 
    202 # width
    203 data = readdlm("width/width_fits.txt")
    204 tau = data[:,2]
    205 s = data[:,3]
    206 plotTau([5e3, 6e3, 7e3, 8e3, 9e3, 1e4],
    207     tau[1:6],
    208     s[1:6],
    209     "Strait width \$W\$ [m]", "width",
    210     x_lim=[4.5e3, 9e3],
    211     y_lim=[1e0, 3e2],
    212     ycenter=2e1,
    213     shadex=[7.2e3, 9e3],
    214     tau2=tau[7:end],
    215     tau2_s=s[7:end],
    216     logx=false)
    217 
    218 # width_monodisperse
    219 data = readdlm("width_monodisperse/width_monodisperse_fits.txt")
    220 tau = data[:,2]
    221 s = data[:,3]
    222 plotTau([6e3, 7e3, 8e3, 9e3, 1e4],
    223     tau[1:5],
    224     s[1:5],
    225     "Strait width \$W\$ [m]", "width_monodisperse",
    226     x_lim=[4e3, 1.1e4],
    227     y_lim=[1e-1, 3e2],
    228     ycenter=3e0,
    229     shadex=[0.1, 0.11],
    230     tau2=tau[6:end],
    231     tau2_s=s[6:end],
    232     logx=false)
    233