Introduction
Introduction Statistics Contact Development Disclaimer Help
tplot.jl - seaice-experiments - sea ice experiments using Granular.jl
git clone git://src.adamsgaard.dk/seaice-experiments
Log
Files
Refs
README
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-fl…
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…
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…
133 #"tau = $(fit.param[1]) h, 95% CI = $(errors[1]…
134 println("$(sim_id) $(fit.param[1]) $(sample_std_dev[1])")
135 #label = @sprintf("\$P_s(t, \\tau = %4.3g\$ h), \$SE_\\tau = %4.…
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…
149 #":C1", linewidth=1)
150 #PyPlot.plot(t_hours, survival_model(t_hours, fit.param + errors…
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_numb…
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()
You are viewing proxied material from mx1.adamsgaard.dk. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.