PySDM_examples.Luettmer_homogeneous_freezing.plot
shared plot functions for homogeneous freezing notebooks
1"""shared plot functions for homogeneous freezing notebooks""" 2 3from matplotlib import pyplot, ticker 4import numpy as np 5import seaborn as sns 6from cycler import cycler 7from PySDM import Formulae 8 9formulae = Formulae(particle_shape_and_density="MixedPhaseSpheres") 10 11ax_title_size = 18 12ax_lab_fsize = 15 13tick_fsize = 15 14T_frz_bins = np.linspace(-40, -34, num=60, endpoint=True) 15T_frz_bins_kelvin = np.linspace(230, 240, num=100, endpoint=True) 16 17 18def cumulative_histogram(data, bins, reverse=False, density=True): 19 hist, bin_edges = np.histogram(data, bins=bins, density=False) 20 21 if reverse: 22 cum_hist = np.cumsum(hist[::-1])[::-1] 23 cum_hist_0 = cum_hist[0] 24 else: 25 cum_hist = np.cumsum(hist) 26 cum_hist_0 = cum_hist[-1] 27 28 if density: 29 cum_hist = cum_hist / cum_hist_0 30 31 bin_centers = 0.5 * (bin_edges[:-1] + bin_edges[1:]) 32 return cum_hist, bin_centers 33 34 35def plot_thermodynamics_and_bulk( 36 simulation, 37 title_add=None, 38 show_conc=False, 39 show_jhom=True, 40 show_tf=True, 41 t_lim=None, 42): 43 if title_add is None: 44 title_add = ["", "", "", ""] 45 plot_daw = False 46 output = simulation["ensemble_member_outputs"][0] 47 time = output["t"] 48 T = np.asarray(output["T"]) 49 RH = np.asarray(output["RH"]) / 100 50 RHi = np.asarray(output["RH_ice"]) / 100 51 qc = np.asarray(output["LWC"]) 52 qi = np.asarray(output["IWC"]) 53 qv = np.asarray(output["qv"]) 54 T_frz = np.asarray(output["T_frz"]) 55 if show_conc: 56 nc, ni = np.asarray(output["ns"]), np.asarray(output["ni"]) 57 58 if t_lim is None: 59 t_lim = np.amax(time) 60 61 first_ice_idx = np.where(qi > 1e-7)[0][0] 62 first_ice_time = time[first_ice_idx] 63 first_T_frz = T[first_ice_idx] 64 65 if not show_tf: 66 rc, ri = np.asarray(output["rs"]), np.asarray(output["ri"]) 67 68 if show_jhom: 69 svp = Formulae( 70 saturation_vapour_pressure="FlatauWalkoCotton" 71 ).saturation_vapour_pressure 72 a_w_ice = svp.pvs_ice(T) / svp.pvs_water(T) 73 d_a_w_ice = (RHi - 1) * a_w_ice 74 75 j_hom_rate = Formulae( 76 homogeneous_ice_nucleation_rate="KoopMurray2016" 77 ).homogeneous_ice_nucleation_rate 78 koop_murray_2016 = j_hom_rate.j_hom(T, d_a_w_ice) 79 j_hom_rate = Formulae( 80 homogeneous_ice_nucleation_rate="Koop_Correction" 81 ).homogeneous_ice_nucleation_rate 82 spichtinger_2023 = j_hom_rate.j_hom(T, d_a_w_ice) 83 abs_diff_j_hom = (koop_murray_2016 - spichtinger_2023) / koop_murray_2016 84 else: 85 radius = np.asarray(output["radius"]) 86 multiplicity = np.asarray(output["multiplicity"]) 87 88 _, axs = pyplot.subplots(1, 4, figsize=(20, 5), constrained_layout=True) 89 90 # Temperture profile 91 iax = 0 92 ax = axs[iax] 93 ax.plot(time, RH, color="red", linestyle="dashdot", label=r"$S_\text{w}$") 94 ax.plot(time, RHi, color="blue", linestyle="--", label=r"$S_\text{i}$") 95 ax.set_ylabel("saturation ratio", fontsize=ax_lab_fsize) 96 ax.legend(loc="center left", fontsize=ax_lab_fsize) 97 ax.set_xlim(time[0], t_lim) 98 ax.tick_params(labelsize=tick_fsize) 99 ax.set_title(title_add[iax] + r"ambient thermodynamics", fontsize=ax_lab_fsize) 100 ax.grid(visible=True) 101 ax.axvline(x=first_ice_time, color="black", linestyle=":") 102 103 twin = ax.twinx() 104 twin.plot(time, T, color="black", linestyle="-", label="T") 105 twin.set_xlabel("time [s]", fontsize=ax_lab_fsize) 106 twin.set_ylabel("temperature [K]", fontsize=ax_lab_fsize) 107 twin.legend(loc="upper left", fontsize=ax_lab_fsize) 108 twin.tick_params(labelsize=tick_fsize) 109 110 # mixing ratio and number concentration 111 iax = 1 112 ax = axs[iax] 113 ax.plot(time, qc, color="red", linestyle="dashdot", label=r"$q_\text{w}$") 114 ax.plot(time, qi, color="blue", linestyle="--", label=r"$q_\text{i}$") 115 ax.plot(time, qv, color="black", linestyle="-", label=r"$q_\text{v}$") 116 ax.set_yscale("log") 117 ax.set_ylim(1e-5, 1e-2) 118 ax.set_xlabel("time [s]", fontsize=ax_lab_fsize) 119 ax.set_ylabel(r"mixing ratio [$\mathrm{kg \, kg^{-1}}$]", fontsize=ax_lab_fsize) 120 ax.legend(fontsize=ax_lab_fsize) 121 ax.tick_params(labelsize=tick_fsize) 122 ax.set_xlim(time[0], t_lim) 123 ax.grid(visible=True) 124 ax.axvline(x=first_ice_time, color="black", linestyle=":") 125 ax.set_title(title_add[iax] + r"bulk quantities", fontsize=ax_lab_fsize) 126 if show_conc: 127 twin = ax.twinx() 128 twin.plot(time, nc, color="red", linestyle="densly dashdot", label="water") 129 twin.plot(time, ni, color="blue", linestyle="densly dashed", label="ice") 130 twin.set_yscale("log") 131 twin.set_xlabel("time [s]", fontsize=ax_lab_fsize) 132 twin.set_ylabel( 133 r"number concentration [$\mathrm{kg^{-1}}$]", fontsize=ax_lab_fsize 134 ) 135 twin.tick_params(labelsize=tick_fsize) 136 137 iax = 2 138 ax = axs[iax] 139 if show_tf: 140 ax.hist( 141 T_frz, 142 bins=T_frz_bins_kelvin, 143 density=True, 144 cumulative=-1, 145 alpha=1.0, 146 histtype="step", 147 linewidth=1.5, 148 ) 149 ax.set_xlim(left=234, right=239) 150 ax.set_xlabel("freezing temperature [K]", fontsize=ax_lab_fsize) 151 ax.set_ylabel("frozen fraction", fontsize=ax_lab_fsize) 152 ax.tick_params(labelsize=tick_fsize) 153 ax.grid(visible=True) 154 ax.axvline(x=first_T_frz, color="black", linestyle=":") 155 ax.set_title(title_add[iax] + r"$T_{frz}$ histogram", fontsize=ax_lab_fsize) 156 else: 157 ax.plot(time, rc * 1e6, color="red", linestyle="dashdot", label="water") 158 ax.plot(time, ri * 1e6, color="blue", linestyle="--", label="ice") 159 ax.set_yscale("log") 160 ax.set_ylim(1e-2, 1e2) 161 ax.set_xlabel("time [s]", fontsize=ax_lab_fsize) 162 ax.set_ylabel("mean radius [µm]", fontsize=ax_lab_fsize) 163 ax.legend(fontsize=ax_lab_fsize) 164 ax.set_xlim(time[0], t_lim) 165 ax.tick_params(labelsize=tick_fsize) 166 ax.grid(visible=True) 167 ax.axvline(x=first_ice_time, color="black", linestyle=":") 168 ax.set_title(title_add[iax] + r" mean radius", fontsize=ax_lab_fsize) 169 170 # Water activity difference profile 171 iax = 3 172 ax = axs[iax] 173 if show_jhom: 174 lin_s_SP2023 = "--" 175 lin_s_KM2016 = "-" 176 if simulation["settings"]["hom_freezing"] == "Spichtinger2023": 177 lin_s_SP2023 = "-" 178 lin_s_KM2016 = "--" 179 180 ax.plot( 181 time, 182 koop_murray_2016, 183 color="blue", 184 linestyle=lin_s_KM2016, 185 label="JHOM-T", 186 ) 187 ax.plot( 188 time, 189 spichtinger_2023, 190 color="red", 191 linestyle=lin_s_SP2023, 192 label="JHOM-DWA", 193 ) 194 ax.set_ylabel( 195 r"nucleation rate [$\mathrm{m^{-3} \, s^{-1}}$]", fontsize=ax_lab_fsize 196 ) 197 ax.set_ylim(1e-30, 1e30) 198 ax.set_title(title_add[iax] + r"nucleation rates", fontsize=ax_lab_fsize) 199 ax.legend(loc="upper left", fontsize=ax_lab_fsize) 200 ax.set_yscale("log") 201 ax.set_xlim(time[0], t_lim) 202 ax.axvline(x=first_ice_time, color="black", linestyle=":") 203 ax.set_xlabel("time [s]", fontsize=ax_lab_fsize) 204 if plot_daw: 205 twin = ax.twinx() 206 twin.plot( 207 time, 208 d_a_w_ice, 209 color="gray", 210 linestyle="dashdot", 211 label=r"$\Delta a_{w}$", 212 ) 213 twin.set_ylim(0.2, 0.35) 214 twin.set_ylabel("water activity difference", fontsize=ax_lab_fsize) 215 twin.tick_params(labelsize=tick_fsize) 216 twin.legend(loc="lower left", fontsize=ax_lab_fsize) 217 else: 218 twin = ax.twinx() 219 twin.plot( 220 time, 221 abs_diff_j_hom, 222 label=r"$\Delta J_{\mathrm{hom}}$", 223 color="gray", 224 linestyle="dashdot", 225 ) 226 twin.set_ylim(-10, 10) 227 twin.set_ylabel("relative error", fontsize=ax_lab_fsize) 228 twin.tick_params(labelsize=tick_fsize) 229 twin.legend(loc="lower left", fontsize=ax_lab_fsize) 230 else: 231 ax.scatter(radius * 1e6, multiplicity) 232 ax.set_yscale("log") 233 ax.set_xscale("log") 234 ax.set_xlim(1e-3, 5e-0) 235 ax.set_xlabel("initial radius [µm]", fontsize=ax_lab_fsize) 236 ax.set_ylabel("multiplicity", fontsize=ax_lab_fsize) 237 ax.set_title(title_add[iax] + r"CCN size distribution", fontsize=ax_lab_fsize) 238 239 ax.tick_params(labelsize=tick_fsize) 240 ax.grid(visible=True) 241 242 243def plot_freezing_temperatures_histogram(ax, simulation, plot_rhi=False): 244 245 number_of_ensemble_runs = simulation["settings"]["number_of_ensemble_runs"] 246 247 for i in range(number_of_ensemble_runs): 248 output = simulation["ensemble_member_outputs"][i] 249 if plot_rhi: 250 var = np.asarray(output["RHi_frz"]) 251 bins = np.linspace(1, 1.6, num=60, endpoint=True) 252 cumulative = 1 253 else: 254 var = np.asarray(output["T_frz"]) 255 bins = T_frz_bins_kelvin 256 cumulative = -1 257 title = "Nucleation rate=" + simulation["settings"]["hom_freezing"] 258 259 ax.hist( 260 var, 261 bins=bins, 262 density=True, 263 cumulative=cumulative, 264 alpha=1.0, 265 histtype="step", 266 linewidth=1.5, 267 ) 268 269 if plot_rhi: 270 ax.set_xlim(left=1.0, right=1.6) 271 ax.set_xlabel("freezing supersaturation wrt ice", fontsize=ax_lab_fsize) 272 else: 273 ax.set_xlim(left=234, right=239) 274 ax.axvline(x=235, color="k", linestyle="--") 275 ax.set_xlabel("freezing temperature [K]", fontsize=ax_lab_fsize) 276 ax.set_title(title, fontsize=ax_lab_fsize) 277 ax.set_ylabel("frozen fraction", fontsize=ax_lab_fsize) 278 ax.tick_params(labelsize=tick_fsize) 279 return ax 280 281 282def plot_freezing_temperatures_histogram_allinone( 283 ax, simulations, title, lloc="upper right" 284): 285 286 colors = ["black", "blue", "red"] 287 linestyles = ["-", "--", ":"] 288 289 for k, simulation in enumerate(simulations): 290 291 number_of_ensemble_runs = simulation["settings"]["number_of_ensemble_runs"] 292 n_sd = simulation["settings"]["n_sd"] 293 histogram_list = np.zeros((number_of_ensemble_runs, len(T_frz_bins_kelvin) - 1)) 294 for i in range(number_of_ensemble_runs): 295 output = simulation["ensemble_member_outputs"][i] 296 T_frz = np.asarray(output["T_frz"]) 297 298 hist, T_frz_bins_center = cumulative_histogram( 299 T_frz, T_frz_bins_kelvin, reverse=True 300 ) 301 histogram_list[i, :] = hist 302 303 max_line = np.max(histogram_list, axis=0) 304 mean_line = np.mean(histogram_list, axis=0) 305 min_line = np.min(histogram_list, axis=0) 306 307 ax.plot( 308 T_frz_bins_center, 309 mean_line, 310 color=colors[k], 311 linestyle=linestyles[k], 312 label=r"$n_\text{sd}$: " + f"{int(n_sd):5.0f}", 313 ) 314 ax.fill_between( 315 T_frz_bins_center, min_line, max_line, color=colors[k], alpha=0.2 316 ) 317 ax.set_xlim(left=234.5, right=239) 318 ax.set_title(title, fontsize=ax_lab_fsize) 319 ax.set_xlabel("freezing temperature [K]", fontsize=ax_lab_fsize) 320 ax.set_ylabel("frozen fraction", fontsize=ax_lab_fsize) 321 ax.tick_params(labelsize=tick_fsize) 322 ax.xaxis.set_major_locator(ticker.MultipleLocator(1)) 323 ax.xaxis.set_minor_locator(ticker.MultipleLocator(0.5)) 324 ax.grid(visible=True) 325 ax.legend(loc=lloc, fontsize=ax_lab_fsize) 326 return ax 327 328 329def plot_freezing_temperatures_2d_histogram_seaborn( 330 ensemble_simulations, 331 hom_freezing_type, 332 title="", 333 height=4, 334 width=5, 335): 336 337 sns.set_theme(style="ticks") 338 339 second_axis = True 340 y_log = True 341 342 ens_variable_name = ensemble_simulations["ens_variable_name"] 343 simulations = ensemble_simulations[hom_freezing_type] 344 345 ens_variable = [] 346 ens_variable_sec = [] 347 T_frz_hist = [] 348 multiplicity_hist = [] 349 350 for simulation in simulations: 351 if ens_variable_name == "sig": 352 ens_variable_value = simulation["settings"]["sigma_droplet_distribution"] 353 else: 354 ens_variable_value = simulation["settings"][ens_variable_name] 355 356 output = simulation["ensemble_member_outputs"][0] 357 qi = np.asarray(output["IWC"]) 358 T_frz = np.asarray(output["T_frz"]) 359 multiplicity = np.asarray(output["multiplicity"]) 360 first_ice_idx = np.where(qi > 1e-7)[0][0] 361 362 if ens_variable_name == "n_ccn": 363 rs = np.asarray(output["rs"]) 364 ens_variable_sec_value = rs[first_ice_idx - 1] 365 else: 366 time = np.asarray(output["t"]) 367 T = np.asarray(output["T"]) 368 ens_variable_sec_value = abs(np.mean(np.diff(T) / np.diff(time))) 369 370 T_frz_hist.extend(T_frz) 371 multiplicity_hist.extend(multiplicity) 372 ens_variable.extend(np.full_like(T_frz, ens_variable_value)) 373 ens_variable_sec.append(ens_variable_sec_value) 374 375 ens_variable = np.array(ens_variable) 376 ens_variable_sec = np.array(ens_variable_sec) 377 y_label, y_label_sec = "", "" 378 if ens_variable_name == "w_updraft": 379 y_label = r"vertical updraft [$\mathrm{m \, s^{-1}}$]" 380 y_label_sec = r"cooling rate [$\mathrm{mK \, s^{-1}}$]" 381 ens_variable_sec = ens_variable_sec * 1e3 382 binwidth = 0.25 383 elif ens_variable_name == "n_ccn": 384 y_label = r"ccn concentration [$\mathrm{cm^{-3}}$]" 385 y_label_sec = r"radius [$\mathrm{\mu m}$]" 386 ens_variable = ens_variable / 1.0e6 387 ens_variable_sec = ens_variable_sec * 1.0e6 388 binwidth = 0.25 389 elif ens_variable_name == "sig": 390 y_label = r"$\sigma$" 391 second_axis = False 392 y_log = False 393 binwidth = 0.1 394 395 ens_variable_label = np.unique(np.sort(ens_variable)) 396 397 xlim = (232.5, 240) 398 h = sns.JointGrid( 399 x=T_frz_hist, 400 y=ens_variable, 401 xlim=xlim, 402 ) 403 if y_log: 404 h.ax_joint.set(yscale="log") 405 406 h.plot_joint( 407 sns.histplot, 408 stat="probability", 409 binwidth=binwidth, 410 discrete=(False, False), 411 weights=multiplicity_hist, 412 pmax=0.8, 413 ) 414 415 h.plot_marginals( 416 sns.histplot, 417 element="step", 418 ) 419 h.set_axis_labels("freezing temperature [K]", y_label, fontsize=ax_lab_fsize) 420 h.ax_joint.set_title( 421 title, 422 pad=50, 423 fontsize=ax_title_size, 424 ) 425 h.ax_joint.tick_params(labelsize=tick_fsize) 426 h.ax_joint.grid(True, axis="x", which="both", linestyle="-", linewidth=0.6) 427 h.ax_marg_x.tick_params(labelsize=tick_fsize) 428 h.ax_marg_y.tick_params(labelsize=tick_fsize) 429 430 h.ax_joint.xaxis.set_major_locator(ticker.MultipleLocator(2, offset=1)) 431 h.ax_joint.xaxis.set_minor_locator(ticker.MultipleLocator(1, offset=1)) 432 h.ax_marg_y.remove() 433 434 if second_axis: 435 436 ax2 = h.ax_joint.secondary_yaxis("right", functions=(lambda y: y, lambda y: y)) 437 ax2.set_yticks(ens_variable_label) 438 ax2.set_yticklabels([f"{v:.1f}" for v in ens_variable_sec]) 439 ax2.minorticks_off() 440 ax2.tick_params( 441 axis="y", 442 which="both", 443 direction="out", 444 length=h.ax_joint.yaxis.get_ticklines()[0].get_markersize(), 445 width=h.ax_joint.yaxis.get_ticklines()[0].get_markeredgewidth(), 446 labelsize=tick_fsize, 447 ) 448 ax2.set_ylabel(y_label_sec, fontsize=ax_lab_fsize) 449 450 h.fig.set_size_inches(width, height) 451 452 453def plot_ensemble_bulk( 454 ax, ensemble_simulations, var_name, title_add="" 455): # pylint: disable=too-many-nested-blocks 456 457 colors = ["blue", "red", "cyan"] 458 linestyles = ["-", "--", ":"] 459 pyplot.rcParams["axes.prop_cycle"] = cycler(color=colors) + cycler( 460 linestyle=linestyles 461 ) 462 463 for ensemble_simulation in ensemble_simulations: 464 ens_var = np.asarray(ensemble_simulation["ens_variable"]) 465 ens_var_name = ensemble_simulation["ens_variable_name"] 466 hom_freezing_types = ensemble_simulation["hom_freezing_types"] 467 hom_freezing_labels = ensemble_simulation["hom_freezing_labels"] 468 len_ens_var = len(ens_var) 469 470 if ens_var_name == "n_ccn": 471 ens_var_scale = 1.0 / 1e6 472 else: 473 ens_var_scale = 1.0 474 if ens_var_name == "sig": 475 ens_var_name = "sigma_droplet_distribution" 476 477 for j, hom_freezing_type in enumerate(hom_freezing_types): 478 simulations = ensemble_simulation[hom_freezing_type] 479 number_of_ensemble_runs = simulations[0]["settings"][ 480 "number_of_ensemble_runs" 481 ] 482 var = np.zeros((len_ens_var, number_of_ensemble_runs)) 483 for i in range(len_ens_var): 484 for simulation in simulations: 485 if simulation["settings"][ens_var_name] == ens_var[i]: 486 for h in range(number_of_ensemble_runs): 487 output = simulation["ensemble_member_outputs"][h] 488 if var_name == "freezing_fraction": 489 ni = np.asarray(output["ni"])[-1] 490 nc = np.asarray(output["ns"])[0] 491 var[i, h] = (1 - (nc - ni) / nc) * 100 492 else: 493 var[i, h] = np.asarray(output[var_name])[-1] 494 495 if number_of_ensemble_runs > 1: 496 ax.fill_betweenx( 497 ens_var * ens_var_scale, 498 np.min(var, axis=1), 499 np.max(var, axis=1), 500 alpha=0.2, 501 ) 502 else: 503 ax.scatter( 504 var[:, 0], 505 ens_var * ens_var_scale, 506 ) 507 ax.plot( 508 var[:, 0], 509 ens_var * ens_var_scale, 510 label=hom_freezing_labels[j], 511 ) 512 513 title, x_label, y_label, ens_label = "", "", "", "" 514 if var_name == "ni": 515 ax.set_xscale("log") 516 x_label = r"ice number concentration [$\mathrm{kg^{-1}}$]" 517 title = r"$n_\text{i}$" 518 ax.set_xlim(1e6, 1e10) 519 elif var_name == "IWC": 520 ax.set_xscale("log") 521 x_label = r"mixing ratio [$\mathrm{kg \, kg^{-1}}$]" 522 title = "ice mixing ratio" 523 ax.set_xlim(1e-4, 1e-3) 524 elif var_name == "freezing_fraction": 525 title = r"$n_\text{frz}$" 526 x_label = r"frozen fraction [$\mathrm{\%}$]" 527 ax.set_xlim(0, 20) 528 529 if ens_var_name == "n_ccn": 530 ax.set_yscale("log") 531 y_label = r"ccn concentration [$\mathrm{cm^{-3}}$]" 532 ens_label = r"$n_\text{ccn}$ ensemble" 533 elif ens_var_name == "w_updraft": 534 ax.set_yscale("log") 535 y_label = r"vertical updraft [$\mathrm{m \, s^{-1}}$]" 536 ens_label = "w ensemble" 537 elif ens_var_name == "sigma_droplet_distribution": 538 y_label = r"standard deviation DSD" 539 ens_label = r"$\sigma$ ensemble" 540 elif ens_var_name == "n_sd": 541 ax.set_yscale("log") 542 y_label = "number of super-particles" 543 ens_label = r"$n_\text{sd}$ ensemble" 544 ax.set_title(title_add + " " + title + " for " + ens_label, fontsize=ax_lab_fsize) 545 ax.set_xlabel(x_label, fontsize=ax_lab_fsize) 546 ax.set_ylabel(y_label, fontsize=ax_lab_fsize) 547 ax.grid(visible=True) 548 ax.legend(fontsize=ax_lab_fsize) 549 return ax
formulae =
<PySDM.formulae.Formulae object>
ax_title_size =
18
ax_lab_fsize =
15
tick_fsize =
15
T_frz_bins =
array([-40. , -39.89830508, -39.79661017, -39.69491525,
-39.59322034, -39.49152542, -39.38983051, -39.28813559,
-39.18644068, -39.08474576, -38.98305085, -38.88135593,
-38.77966102, -38.6779661 , -38.57627119, -38.47457627,
-38.37288136, -38.27118644, -38.16949153, -38.06779661,
-37.96610169, -37.86440678, -37.76271186, -37.66101695,
-37.55932203, -37.45762712, -37.3559322 , -37.25423729,
-37.15254237, -37.05084746, -36.94915254, -36.84745763,
-36.74576271, -36.6440678 , -36.54237288, -36.44067797,
-36.33898305, -36.23728814, -36.13559322, -36.03389831,
-35.93220339, -35.83050847, -35.72881356, -35.62711864,
-35.52542373, -35.42372881, -35.3220339 , -35.22033898,
-35.11864407, -35.01694915, -34.91525424, -34.81355932,
-34.71186441, -34.61016949, -34.50847458, -34.40677966,
-34.30508475, -34.20338983, -34.10169492, -34. ])
T_frz_bins_kelvin =
array([230. , 230.1010101 , 230.2020202 , 230.3030303 ,
230.4040404 , 230.50505051, 230.60606061, 230.70707071,
230.80808081, 230.90909091, 231.01010101, 231.11111111,
231.21212121, 231.31313131, 231.41414141, 231.51515152,
231.61616162, 231.71717172, 231.81818182, 231.91919192,
232.02020202, 232.12121212, 232.22222222, 232.32323232,
232.42424242, 232.52525253, 232.62626263, 232.72727273,
232.82828283, 232.92929293, 233.03030303, 233.13131313,
233.23232323, 233.33333333, 233.43434343, 233.53535354,
233.63636364, 233.73737374, 233.83838384, 233.93939394,
234.04040404, 234.14141414, 234.24242424, 234.34343434,
234.44444444, 234.54545455, 234.64646465, 234.74747475,
234.84848485, 234.94949495, 235.05050505, 235.15151515,
235.25252525, 235.35353535, 235.45454545, 235.55555556,
235.65656566, 235.75757576, 235.85858586, 235.95959596,
236.06060606, 236.16161616, 236.26262626, 236.36363636,
236.46464646, 236.56565657, 236.66666667, 236.76767677,
236.86868687, 236.96969697, 237.07070707, 237.17171717,
237.27272727, 237.37373737, 237.47474747, 237.57575758,
237.67676768, 237.77777778, 237.87878788, 237.97979798,
238.08080808, 238.18181818, 238.28282828, 238.38383838,
238.48484848, 238.58585859, 238.68686869, 238.78787879,
238.88888889, 238.98989899, 239.09090909, 239.19191919,
239.29292929, 239.39393939, 239.49494949, 239.5959596 ,
239.6969697 , 239.7979798 , 239.8989899 , 240. ])
def
cumulative_histogram(data, bins, reverse=False, density=True):
19def cumulative_histogram(data, bins, reverse=False, density=True): 20 hist, bin_edges = np.histogram(data, bins=bins, density=False) 21 22 if reverse: 23 cum_hist = np.cumsum(hist[::-1])[::-1] 24 cum_hist_0 = cum_hist[0] 25 else: 26 cum_hist = np.cumsum(hist) 27 cum_hist_0 = cum_hist[-1] 28 29 if density: 30 cum_hist = cum_hist / cum_hist_0 31 32 bin_centers = 0.5 * (bin_edges[:-1] + bin_edges[1:]) 33 return cum_hist, bin_centers
def
plot_thermodynamics_and_bulk( simulation, title_add=None, show_conc=False, show_jhom=True, show_tf=True, t_lim=None):
36def plot_thermodynamics_and_bulk( 37 simulation, 38 title_add=None, 39 show_conc=False, 40 show_jhom=True, 41 show_tf=True, 42 t_lim=None, 43): 44 if title_add is None: 45 title_add = ["", "", "", ""] 46 plot_daw = False 47 output = simulation["ensemble_member_outputs"][0] 48 time = output["t"] 49 T = np.asarray(output["T"]) 50 RH = np.asarray(output["RH"]) / 100 51 RHi = np.asarray(output["RH_ice"]) / 100 52 qc = np.asarray(output["LWC"]) 53 qi = np.asarray(output["IWC"]) 54 qv = np.asarray(output["qv"]) 55 T_frz = np.asarray(output["T_frz"]) 56 if show_conc: 57 nc, ni = np.asarray(output["ns"]), np.asarray(output["ni"]) 58 59 if t_lim is None: 60 t_lim = np.amax(time) 61 62 first_ice_idx = np.where(qi > 1e-7)[0][0] 63 first_ice_time = time[first_ice_idx] 64 first_T_frz = T[first_ice_idx] 65 66 if not show_tf: 67 rc, ri = np.asarray(output["rs"]), np.asarray(output["ri"]) 68 69 if show_jhom: 70 svp = Formulae( 71 saturation_vapour_pressure="FlatauWalkoCotton" 72 ).saturation_vapour_pressure 73 a_w_ice = svp.pvs_ice(T) / svp.pvs_water(T) 74 d_a_w_ice = (RHi - 1) * a_w_ice 75 76 j_hom_rate = Formulae( 77 homogeneous_ice_nucleation_rate="KoopMurray2016" 78 ).homogeneous_ice_nucleation_rate 79 koop_murray_2016 = j_hom_rate.j_hom(T, d_a_w_ice) 80 j_hom_rate = Formulae( 81 homogeneous_ice_nucleation_rate="Koop_Correction" 82 ).homogeneous_ice_nucleation_rate 83 spichtinger_2023 = j_hom_rate.j_hom(T, d_a_w_ice) 84 abs_diff_j_hom = (koop_murray_2016 - spichtinger_2023) / koop_murray_2016 85 else: 86 radius = np.asarray(output["radius"]) 87 multiplicity = np.asarray(output["multiplicity"]) 88 89 _, axs = pyplot.subplots(1, 4, figsize=(20, 5), constrained_layout=True) 90 91 # Temperture profile 92 iax = 0 93 ax = axs[iax] 94 ax.plot(time, RH, color="red", linestyle="dashdot", label=r"$S_\text{w}$") 95 ax.plot(time, RHi, color="blue", linestyle="--", label=r"$S_\text{i}$") 96 ax.set_ylabel("saturation ratio", fontsize=ax_lab_fsize) 97 ax.legend(loc="center left", fontsize=ax_lab_fsize) 98 ax.set_xlim(time[0], t_lim) 99 ax.tick_params(labelsize=tick_fsize) 100 ax.set_title(title_add[iax] + r"ambient thermodynamics", fontsize=ax_lab_fsize) 101 ax.grid(visible=True) 102 ax.axvline(x=first_ice_time, color="black", linestyle=":") 103 104 twin = ax.twinx() 105 twin.plot(time, T, color="black", linestyle="-", label="T") 106 twin.set_xlabel("time [s]", fontsize=ax_lab_fsize) 107 twin.set_ylabel("temperature [K]", fontsize=ax_lab_fsize) 108 twin.legend(loc="upper left", fontsize=ax_lab_fsize) 109 twin.tick_params(labelsize=tick_fsize) 110 111 # mixing ratio and number concentration 112 iax = 1 113 ax = axs[iax] 114 ax.plot(time, qc, color="red", linestyle="dashdot", label=r"$q_\text{w}$") 115 ax.plot(time, qi, color="blue", linestyle="--", label=r"$q_\text{i}$") 116 ax.plot(time, qv, color="black", linestyle="-", label=r"$q_\text{v}$") 117 ax.set_yscale("log") 118 ax.set_ylim(1e-5, 1e-2) 119 ax.set_xlabel("time [s]", fontsize=ax_lab_fsize) 120 ax.set_ylabel(r"mixing ratio [$\mathrm{kg \, kg^{-1}}$]", fontsize=ax_lab_fsize) 121 ax.legend(fontsize=ax_lab_fsize) 122 ax.tick_params(labelsize=tick_fsize) 123 ax.set_xlim(time[0], t_lim) 124 ax.grid(visible=True) 125 ax.axvline(x=first_ice_time, color="black", linestyle=":") 126 ax.set_title(title_add[iax] + r"bulk quantities", fontsize=ax_lab_fsize) 127 if show_conc: 128 twin = ax.twinx() 129 twin.plot(time, nc, color="red", linestyle="densly dashdot", label="water") 130 twin.plot(time, ni, color="blue", linestyle="densly dashed", label="ice") 131 twin.set_yscale("log") 132 twin.set_xlabel("time [s]", fontsize=ax_lab_fsize) 133 twin.set_ylabel( 134 r"number concentration [$\mathrm{kg^{-1}}$]", fontsize=ax_lab_fsize 135 ) 136 twin.tick_params(labelsize=tick_fsize) 137 138 iax = 2 139 ax = axs[iax] 140 if show_tf: 141 ax.hist( 142 T_frz, 143 bins=T_frz_bins_kelvin, 144 density=True, 145 cumulative=-1, 146 alpha=1.0, 147 histtype="step", 148 linewidth=1.5, 149 ) 150 ax.set_xlim(left=234, right=239) 151 ax.set_xlabel("freezing temperature [K]", fontsize=ax_lab_fsize) 152 ax.set_ylabel("frozen fraction", fontsize=ax_lab_fsize) 153 ax.tick_params(labelsize=tick_fsize) 154 ax.grid(visible=True) 155 ax.axvline(x=first_T_frz, color="black", linestyle=":") 156 ax.set_title(title_add[iax] + r"$T_{frz}$ histogram", fontsize=ax_lab_fsize) 157 else: 158 ax.plot(time, rc * 1e6, color="red", linestyle="dashdot", label="water") 159 ax.plot(time, ri * 1e6, color="blue", linestyle="--", label="ice") 160 ax.set_yscale("log") 161 ax.set_ylim(1e-2, 1e2) 162 ax.set_xlabel("time [s]", fontsize=ax_lab_fsize) 163 ax.set_ylabel("mean radius [µm]", fontsize=ax_lab_fsize) 164 ax.legend(fontsize=ax_lab_fsize) 165 ax.set_xlim(time[0], t_lim) 166 ax.tick_params(labelsize=tick_fsize) 167 ax.grid(visible=True) 168 ax.axvline(x=first_ice_time, color="black", linestyle=":") 169 ax.set_title(title_add[iax] + r" mean radius", fontsize=ax_lab_fsize) 170 171 # Water activity difference profile 172 iax = 3 173 ax = axs[iax] 174 if show_jhom: 175 lin_s_SP2023 = "--" 176 lin_s_KM2016 = "-" 177 if simulation["settings"]["hom_freezing"] == "Spichtinger2023": 178 lin_s_SP2023 = "-" 179 lin_s_KM2016 = "--" 180 181 ax.plot( 182 time, 183 koop_murray_2016, 184 color="blue", 185 linestyle=lin_s_KM2016, 186 label="JHOM-T", 187 ) 188 ax.plot( 189 time, 190 spichtinger_2023, 191 color="red", 192 linestyle=lin_s_SP2023, 193 label="JHOM-DWA", 194 ) 195 ax.set_ylabel( 196 r"nucleation rate [$\mathrm{m^{-3} \, s^{-1}}$]", fontsize=ax_lab_fsize 197 ) 198 ax.set_ylim(1e-30, 1e30) 199 ax.set_title(title_add[iax] + r"nucleation rates", fontsize=ax_lab_fsize) 200 ax.legend(loc="upper left", fontsize=ax_lab_fsize) 201 ax.set_yscale("log") 202 ax.set_xlim(time[0], t_lim) 203 ax.axvline(x=first_ice_time, color="black", linestyle=":") 204 ax.set_xlabel("time [s]", fontsize=ax_lab_fsize) 205 if plot_daw: 206 twin = ax.twinx() 207 twin.plot( 208 time, 209 d_a_w_ice, 210 color="gray", 211 linestyle="dashdot", 212 label=r"$\Delta a_{w}$", 213 ) 214 twin.set_ylim(0.2, 0.35) 215 twin.set_ylabel("water activity difference", fontsize=ax_lab_fsize) 216 twin.tick_params(labelsize=tick_fsize) 217 twin.legend(loc="lower left", fontsize=ax_lab_fsize) 218 else: 219 twin = ax.twinx() 220 twin.plot( 221 time, 222 abs_diff_j_hom, 223 label=r"$\Delta J_{\mathrm{hom}}$", 224 color="gray", 225 linestyle="dashdot", 226 ) 227 twin.set_ylim(-10, 10) 228 twin.set_ylabel("relative error", fontsize=ax_lab_fsize) 229 twin.tick_params(labelsize=tick_fsize) 230 twin.legend(loc="lower left", fontsize=ax_lab_fsize) 231 else: 232 ax.scatter(radius * 1e6, multiplicity) 233 ax.set_yscale("log") 234 ax.set_xscale("log") 235 ax.set_xlim(1e-3, 5e-0) 236 ax.set_xlabel("initial radius [µm]", fontsize=ax_lab_fsize) 237 ax.set_ylabel("multiplicity", fontsize=ax_lab_fsize) 238 ax.set_title(title_add[iax] + r"CCN size distribution", fontsize=ax_lab_fsize) 239 240 ax.tick_params(labelsize=tick_fsize) 241 ax.grid(visible=True)
def
plot_freezing_temperatures_histogram(ax, simulation, plot_rhi=False):
244def plot_freezing_temperatures_histogram(ax, simulation, plot_rhi=False): 245 246 number_of_ensemble_runs = simulation["settings"]["number_of_ensemble_runs"] 247 248 for i in range(number_of_ensemble_runs): 249 output = simulation["ensemble_member_outputs"][i] 250 if plot_rhi: 251 var = np.asarray(output["RHi_frz"]) 252 bins = np.linspace(1, 1.6, num=60, endpoint=True) 253 cumulative = 1 254 else: 255 var = np.asarray(output["T_frz"]) 256 bins = T_frz_bins_kelvin 257 cumulative = -1 258 title = "Nucleation rate=" + simulation["settings"]["hom_freezing"] 259 260 ax.hist( 261 var, 262 bins=bins, 263 density=True, 264 cumulative=cumulative, 265 alpha=1.0, 266 histtype="step", 267 linewidth=1.5, 268 ) 269 270 if plot_rhi: 271 ax.set_xlim(left=1.0, right=1.6) 272 ax.set_xlabel("freezing supersaturation wrt ice", fontsize=ax_lab_fsize) 273 else: 274 ax.set_xlim(left=234, right=239) 275 ax.axvline(x=235, color="k", linestyle="--") 276 ax.set_xlabel("freezing temperature [K]", fontsize=ax_lab_fsize) 277 ax.set_title(title, fontsize=ax_lab_fsize) 278 ax.set_ylabel("frozen fraction", fontsize=ax_lab_fsize) 279 ax.tick_params(labelsize=tick_fsize) 280 return ax
def
plot_freezing_temperatures_histogram_allinone(ax, simulations, title, lloc='upper right'):
283def plot_freezing_temperatures_histogram_allinone( 284 ax, simulations, title, lloc="upper right" 285): 286 287 colors = ["black", "blue", "red"] 288 linestyles = ["-", "--", ":"] 289 290 for k, simulation in enumerate(simulations): 291 292 number_of_ensemble_runs = simulation["settings"]["number_of_ensemble_runs"] 293 n_sd = simulation["settings"]["n_sd"] 294 histogram_list = np.zeros((number_of_ensemble_runs, len(T_frz_bins_kelvin) - 1)) 295 for i in range(number_of_ensemble_runs): 296 output = simulation["ensemble_member_outputs"][i] 297 T_frz = np.asarray(output["T_frz"]) 298 299 hist, T_frz_bins_center = cumulative_histogram( 300 T_frz, T_frz_bins_kelvin, reverse=True 301 ) 302 histogram_list[i, :] = hist 303 304 max_line = np.max(histogram_list, axis=0) 305 mean_line = np.mean(histogram_list, axis=0) 306 min_line = np.min(histogram_list, axis=0) 307 308 ax.plot( 309 T_frz_bins_center, 310 mean_line, 311 color=colors[k], 312 linestyle=linestyles[k], 313 label=r"$n_\text{sd}$: " + f"{int(n_sd):5.0f}", 314 ) 315 ax.fill_between( 316 T_frz_bins_center, min_line, max_line, color=colors[k], alpha=0.2 317 ) 318 ax.set_xlim(left=234.5, right=239) 319 ax.set_title(title, fontsize=ax_lab_fsize) 320 ax.set_xlabel("freezing temperature [K]", fontsize=ax_lab_fsize) 321 ax.set_ylabel("frozen fraction", fontsize=ax_lab_fsize) 322 ax.tick_params(labelsize=tick_fsize) 323 ax.xaxis.set_major_locator(ticker.MultipleLocator(1)) 324 ax.xaxis.set_minor_locator(ticker.MultipleLocator(0.5)) 325 ax.grid(visible=True) 326 ax.legend(loc=lloc, fontsize=ax_lab_fsize) 327 return ax
def
plot_freezing_temperatures_2d_histogram_seaborn(ensemble_simulations, hom_freezing_type, title='', height=4, width=5):
330def plot_freezing_temperatures_2d_histogram_seaborn( 331 ensemble_simulations, 332 hom_freezing_type, 333 title="", 334 height=4, 335 width=5, 336): 337 338 sns.set_theme(style="ticks") 339 340 second_axis = True 341 y_log = True 342 343 ens_variable_name = ensemble_simulations["ens_variable_name"] 344 simulations = ensemble_simulations[hom_freezing_type] 345 346 ens_variable = [] 347 ens_variable_sec = [] 348 T_frz_hist = [] 349 multiplicity_hist = [] 350 351 for simulation in simulations: 352 if ens_variable_name == "sig": 353 ens_variable_value = simulation["settings"]["sigma_droplet_distribution"] 354 else: 355 ens_variable_value = simulation["settings"][ens_variable_name] 356 357 output = simulation["ensemble_member_outputs"][0] 358 qi = np.asarray(output["IWC"]) 359 T_frz = np.asarray(output["T_frz"]) 360 multiplicity = np.asarray(output["multiplicity"]) 361 first_ice_idx = np.where(qi > 1e-7)[0][0] 362 363 if ens_variable_name == "n_ccn": 364 rs = np.asarray(output["rs"]) 365 ens_variable_sec_value = rs[first_ice_idx - 1] 366 else: 367 time = np.asarray(output["t"]) 368 T = np.asarray(output["T"]) 369 ens_variable_sec_value = abs(np.mean(np.diff(T) / np.diff(time))) 370 371 T_frz_hist.extend(T_frz) 372 multiplicity_hist.extend(multiplicity) 373 ens_variable.extend(np.full_like(T_frz, ens_variable_value)) 374 ens_variable_sec.append(ens_variable_sec_value) 375 376 ens_variable = np.array(ens_variable) 377 ens_variable_sec = np.array(ens_variable_sec) 378 y_label, y_label_sec = "", "" 379 if ens_variable_name == "w_updraft": 380 y_label = r"vertical updraft [$\mathrm{m \, s^{-1}}$]" 381 y_label_sec = r"cooling rate [$\mathrm{mK \, s^{-1}}$]" 382 ens_variable_sec = ens_variable_sec * 1e3 383 binwidth = 0.25 384 elif ens_variable_name == "n_ccn": 385 y_label = r"ccn concentration [$\mathrm{cm^{-3}}$]" 386 y_label_sec = r"radius [$\mathrm{\mu m}$]" 387 ens_variable = ens_variable / 1.0e6 388 ens_variable_sec = ens_variable_sec * 1.0e6 389 binwidth = 0.25 390 elif ens_variable_name == "sig": 391 y_label = r"$\sigma$" 392 second_axis = False 393 y_log = False 394 binwidth = 0.1 395 396 ens_variable_label = np.unique(np.sort(ens_variable)) 397 398 xlim = (232.5, 240) 399 h = sns.JointGrid( 400 x=T_frz_hist, 401 y=ens_variable, 402 xlim=xlim, 403 ) 404 if y_log: 405 h.ax_joint.set(yscale="log") 406 407 h.plot_joint( 408 sns.histplot, 409 stat="probability", 410 binwidth=binwidth, 411 discrete=(False, False), 412 weights=multiplicity_hist, 413 pmax=0.8, 414 ) 415 416 h.plot_marginals( 417 sns.histplot, 418 element="step", 419 ) 420 h.set_axis_labels("freezing temperature [K]", y_label, fontsize=ax_lab_fsize) 421 h.ax_joint.set_title( 422 title, 423 pad=50, 424 fontsize=ax_title_size, 425 ) 426 h.ax_joint.tick_params(labelsize=tick_fsize) 427 h.ax_joint.grid(True, axis="x", which="both", linestyle="-", linewidth=0.6) 428 h.ax_marg_x.tick_params(labelsize=tick_fsize) 429 h.ax_marg_y.tick_params(labelsize=tick_fsize) 430 431 h.ax_joint.xaxis.set_major_locator(ticker.MultipleLocator(2, offset=1)) 432 h.ax_joint.xaxis.set_minor_locator(ticker.MultipleLocator(1, offset=1)) 433 h.ax_marg_y.remove() 434 435 if second_axis: 436 437 ax2 = h.ax_joint.secondary_yaxis("right", functions=(lambda y: y, lambda y: y)) 438 ax2.set_yticks(ens_variable_label) 439 ax2.set_yticklabels([f"{v:.1f}" for v in ens_variable_sec]) 440 ax2.minorticks_off() 441 ax2.tick_params( 442 axis="y", 443 which="both", 444 direction="out", 445 length=h.ax_joint.yaxis.get_ticklines()[0].get_markersize(), 446 width=h.ax_joint.yaxis.get_ticklines()[0].get_markeredgewidth(), 447 labelsize=tick_fsize, 448 ) 449 ax2.set_ylabel(y_label_sec, fontsize=ax_lab_fsize) 450 451 h.fig.set_size_inches(width, height)
def
plot_ensemble_bulk(ax, ensemble_simulations, var_name, title_add=''):
454def plot_ensemble_bulk( 455 ax, ensemble_simulations, var_name, title_add="" 456): # pylint: disable=too-many-nested-blocks 457 458 colors = ["blue", "red", "cyan"] 459 linestyles = ["-", "--", ":"] 460 pyplot.rcParams["axes.prop_cycle"] = cycler(color=colors) + cycler( 461 linestyle=linestyles 462 ) 463 464 for ensemble_simulation in ensemble_simulations: 465 ens_var = np.asarray(ensemble_simulation["ens_variable"]) 466 ens_var_name = ensemble_simulation["ens_variable_name"] 467 hom_freezing_types = ensemble_simulation["hom_freezing_types"] 468 hom_freezing_labels = ensemble_simulation["hom_freezing_labels"] 469 len_ens_var = len(ens_var) 470 471 if ens_var_name == "n_ccn": 472 ens_var_scale = 1.0 / 1e6 473 else: 474 ens_var_scale = 1.0 475 if ens_var_name == "sig": 476 ens_var_name = "sigma_droplet_distribution" 477 478 for j, hom_freezing_type in enumerate(hom_freezing_types): 479 simulations = ensemble_simulation[hom_freezing_type] 480 number_of_ensemble_runs = simulations[0]["settings"][ 481 "number_of_ensemble_runs" 482 ] 483 var = np.zeros((len_ens_var, number_of_ensemble_runs)) 484 for i in range(len_ens_var): 485 for simulation in simulations: 486 if simulation["settings"][ens_var_name] == ens_var[i]: 487 for h in range(number_of_ensemble_runs): 488 output = simulation["ensemble_member_outputs"][h] 489 if var_name == "freezing_fraction": 490 ni = np.asarray(output["ni"])[-1] 491 nc = np.asarray(output["ns"])[0] 492 var[i, h] = (1 - (nc - ni) / nc) * 100 493 else: 494 var[i, h] = np.asarray(output[var_name])[-1] 495 496 if number_of_ensemble_runs > 1: 497 ax.fill_betweenx( 498 ens_var * ens_var_scale, 499 np.min(var, axis=1), 500 np.max(var, axis=1), 501 alpha=0.2, 502 ) 503 else: 504 ax.scatter( 505 var[:, 0], 506 ens_var * ens_var_scale, 507 ) 508 ax.plot( 509 var[:, 0], 510 ens_var * ens_var_scale, 511 label=hom_freezing_labels[j], 512 ) 513 514 title, x_label, y_label, ens_label = "", "", "", "" 515 if var_name == "ni": 516 ax.set_xscale("log") 517 x_label = r"ice number concentration [$\mathrm{kg^{-1}}$]" 518 title = r"$n_\text{i}$" 519 ax.set_xlim(1e6, 1e10) 520 elif var_name == "IWC": 521 ax.set_xscale("log") 522 x_label = r"mixing ratio [$\mathrm{kg \, kg^{-1}}$]" 523 title = "ice mixing ratio" 524 ax.set_xlim(1e-4, 1e-3) 525 elif var_name == "freezing_fraction": 526 title = r"$n_\text{frz}$" 527 x_label = r"frozen fraction [$\mathrm{\%}$]" 528 ax.set_xlim(0, 20) 529 530 if ens_var_name == "n_ccn": 531 ax.set_yscale("log") 532 y_label = r"ccn concentration [$\mathrm{cm^{-3}}$]" 533 ens_label = r"$n_\text{ccn}$ ensemble" 534 elif ens_var_name == "w_updraft": 535 ax.set_yscale("log") 536 y_label = r"vertical updraft [$\mathrm{m \, s^{-1}}$]" 537 ens_label = "w ensemble" 538 elif ens_var_name == "sigma_droplet_distribution": 539 y_label = r"standard deviation DSD" 540 ens_label = r"$\sigma$ ensemble" 541 elif ens_var_name == "n_sd": 542 ax.set_yscale("log") 543 y_label = "number of super-particles" 544 ens_label = r"$n_\text{sd}$ ensemble" 545 ax.set_title(title_add + " " + title + " for " + ens_label, fontsize=ax_lab_fsize) 546 ax.set_xlabel(x_label, fontsize=ax_lab_fsize) 547 ax.set_ylabel(y_label, fontsize=ax_lab_fsize) 548 ax.grid(visible=True) 549 ax.legend(fontsize=ax_lab_fsize) 550 return ax