tplot.jl - seaice-experiments - sea ice experiments using Granular.jl
(HTM) git clone git://src.adamsgaard.dk/seaice-experiments
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
tplot.jl (6853B)
---
1 #!/usr/bin/env julia
2 import ArgParse
3 ENV["MPLBACKEND"] = "Agg"
4 import PyPlot
5 import LsqFit
6
7 verbose = false
8
9 function parse_command_line()
10 s = ArgParse.ArgParseSettings()
11 ArgParse.@add_arg_table s begin
12 "--nruns"
13 help = "number of runs in ensemble"
14 arg_type = Int
15 default = 1
16 "id"
17 help = "simulation id"
18 required = true
19 end
20 return ArgParse.parse_args(s)
21 end
22
23 function report_args(parsed_args)
24 println("Parsed args:")
25 for (arg,val) in parsed_args
26 println(" $arg => $val")
27 end
28 end
29
30 function main()
31 parsed_args = parse_command_line()
32 #report_args(parsed_args)
33
34 sim_id = parsed_args["id"]
35 #info(sim_id)
36 nruns = parsed_args["nruns"]
37 t = Vector
38 #t_jam = ones(nruns)*parsed_args["total_hours"]*60.*60.*2.
39 t_jam = ones(nruns)*1e16
40 nskipped = 0
41 iteration_skipped = zeros(Bool, nruns)
42 avg_coordination_number = Vector[]
43
44 PyPlot.figure(figsize=(5,4))
45
46 for i=1:nruns
47 seed = i
48 filename = parsed_args["id"] * "-seed" * string(seed) * "-ice-flux.txt"
49 #info(filename)
50 if !isfile(filename)
51 nskipped += 1
52 iteration_skipped[i] = true
53 push!(avg_coordination_number, [])
54 continue
55 end
56 indata = readdlm(filename)
57 #t, ice_flux = indata[1,:], indata[2,:]
58 t, ice_flux, avg_coordination_no =
59 indata[1,10:end], indata[2,10:end], indata[3,10:end] # skip first 10 vals
60 t -= t[1]
61 push!(avg_coordination_number, avg_coordination_no)
62 PyPlot.plot(t/(60.*60.), ice_flux)
63
64 time_elapsed_while_jammed = 0.
65 for it=(length(t) - 1):-1:1
66 if ice_flux[it] ≈ ice_flux[end]
67 time_elapsed_while_jammed = t[end] - t[it]
68 end
69 end
70
71 if time_elapsed_while_jammed > 60.*60.
72 t_jam[i] = t[end] - time_elapsed_while_jammed
73 #info("simulation $i jammed at t = $(t_jam[i]/(60.*60.)) h")
74 end
75
76 end
77 #PyPlot.title(parsed_args["id"])
78 PyPlot.xlabel("Time [h]")
79 PyPlot.ylabel("Cumulative ice throughput [kg]")
80 PyPlot.tight_layout()
81 PyPlot.title("(a)")
82
83 PyPlot.savefig(parsed_args["id"] * ".png")
84 PyPlot.savefig(parsed_args["id"])
85
86 PyPlot.clf()
87 jam_fraction = zeros(length(t))
88 ensemble_avg_coordination_number = zeros(length(t))
89 for it=1:length(t)
90 for i=1:nruns
91 if iteration_skipped[i]
92 continue
93 end
94 if t_jam[i] <= t[it]
95 jam_fraction[it] += 1./float(nruns - nskipped)
96 end
97 if it > length(avg_coordination_number[i])
98 continue
99 end
100 ensemble_avg_coordination_number[it] +=
101 avg_coordination_number[i][it]/float(nruns - nskipped)
102 end
103 end
104 survived_fraction = 1. - jam_fraction
105 t_hours = t/(60.*60.)
106
107 # fit function of exponential decay to survived_fraction
108 survival_model(t_hours, tau) = exp.(-t_hours/tau[1])
109 tau0 = [2.] # initial guess of tau [hours]
110 I = find(jam_fraction) # find first where survived_fraction > 0.
111 first_idx = 1
112 if length(I) > 0
113 first_idx = I[1]
114 end
115 if t_hours[first_idx] > 6.
116 first_idx = 1
117 end
118 t_hours_trim = linspace(0., t_hours[end],
119 length(survived_fraction[first_idx:end]))
120 fit = LsqFit.curve_fit(survival_model,
121 t_hours_trim,
122 survived_fraction[first_idx:end],
123 tau0)
124
125 covar = LsqFit.estimate_covar(fit)
126 std_error = sqrt.(abs.(diag(covar)))
127 sample_std_dev = std_error .* sqrt(length(t_hours_trim))
128 errors = LsqFit.estimate_errors(fit, 0.95)
129
130 if fit.converged
131 #print_with_color(:green,
132 #"tau = $(fit.param[1]) h, 95% s = $(sample_std_dev[1]) h\n")
133 #"tau = $(fit.param[1]) h, 95% CI = $(errors[1]) h\n")
134 println("$(sim_id) $(fit.param[1]) $(sample_std_dev[1])")
135 #label = @sprintf("\$P_s(t, \\tau = %4.3g\$ h), \$SE_\\tau = %4.3g\$ h",
136 #fit.param[1], std_error[1])
137 label = @sprintf("\$P_s(t, T = %4.3g\$ h), \$s_T = %4.3g\$ h",
138 fit.param[1], sample_std_dev[1])
139
140 # 95% CI range for normal distribution
141 #lower_CI = survival_model(t_hours, fit.param - 1.96*std_error[1])
142 #upper_CI = survival_model(t_hours, fit.param + 1.96*std_error[1])
143
144 # +/- one sample standard deviation
145 lower_SD = survival_model(t_hours, fit.param - sample_std_dev[1])
146 upper_SD = survival_model(t_hours, fit.param + sample_std_dev[1])
147
148 #PyPlot.plot(t_hours, survival_model(t_hours, fit.param - errors[1]),
149 #":C1", linewidth=1)
150 #PyPlot.plot(t_hours, survival_model(t_hours, fit.param + errors[1]),
151 #":C1", linewidth=1)
152
153 PyPlot.fill_between(t_hours + t_hours[first_idx],
154 lower_SD, upper_SD, facecolor="C1",
155 interpolate=true, alpha=0.33)
156
157 PyPlot.plot(t_hours + t_hours[first_idx],
158 survival_model(t_hours, fit.param),
159 "--C1", linewidth=2, label=label)
160
161 else
162 warn(parsed_args["id"] * ": LsqFit.curve_fit did not converge")
163 end
164
165 # Plot DEM
166 PyPlot.plot(t_hours, survived_fraction,
167 "-C0", linewidth=2,
168 label="DEM (\$n = $(nruns - nskipped) \$)")
169
170 if fit.param[1] > 30.
171 PyPlot.legend(loc="lower right", fancybox="true")
172 else
173 PyPlot.legend(loc="upper right", fancybox="true")
174 end
175
176 #PyPlot.title(parsed_args["id"] * ", N = " * string(nruns - nskipped))
177 PyPlot.title("(b)")
178 PyPlot.xlabel("Time [h]")
179 PyPlot.ylabel("Probability of survival \$P_s\$ [-]")
180 PyPlot.xlim([minimum(t_hours), maximum(t_hours)])
181 PyPlot.ylim([-0.02, 1.02])
182 #PyPlot.ylim([0., 1.])
183
184 PyPlot.tight_layout()
185 PyPlot.savefig(parsed_args["id"] * "-survived_fraction.pdf")
186 PyPlot.savefig(parsed_args["id"] * "-survived_fraction.pdf.png")
187 writedlm(parsed_args["id"] * "-survived_fraction.txt",
188 [t_hours, survived_fraction, ensemble_avg_coordination_number])
189
190 # Plot ensemble-mean coordination number
191 PyPlot.clf()
192 PyPlot.plot(t_hours, ensemble_avg_coordination_number,
193 "-C0", linewidth=2)
194 PyPlot.xlabel("Time [h]")
195 PyPlot.ylabel("Avg. coordination number \$Z\$ [-]")
196 PyPlot.xlim([minimum(t_hours), maximum(t_hours)])
197 PyPlot.ylim([-0.02, 5.02])
198 PyPlot.tight_layout()
199 PyPlot.savefig(parsed_args["id"] * "-coordination_number.pdf")
200 PyPlot.savefig(parsed_args["id"] * "-coordination_number.pdf.png")
201
202 end
203
204 main()