pydfnWorks
python wrapper for dfnWorks
plot_fracture_radii.py
Go to the documentation of this file.
1 """
2  :filename: plot_fracture_radii.py
3  :synopsis: Make plots of fracture radii for the entire network and by family. Plots of the expected analytic distributions are included in the plots for each family
4  :version: 1.0
5  :maintainer: Jeffrey Hyman
6  :moduleauthor: Jeffrey Hyman <jhyman@lanl.gov>
7 """
8 
9 import numpy as np
10 import matplotlib.pylab as plt
11 
12 from pydfnworks.dfnGen.generation.output_report.distributions import create_ecdf, tpl, lognormal, exponential
13 
14 
15 def plot_fracture_radii(params, families, fractures, num_bins=20):
16  """ Creates histogram plots of fracture centers in the domain. First all fractures are plotted, then by family.
17 
18  Parameters
19  ------------
20  params : dictionary
21  General dictionary of output analysis code. Contains information on number of families, number of fractures, and colors.
22 
23  families: list of fracture family dictionaries
24  Created by get_family_information
25 
26  fractures: list of fracture dictionaries
27  Created by get_fracture_information
28 
29  Returns
30  --------
31  None
32 
33  Notes
34  -------
35  PDF files are dumped into family/figures. There is one figure for all fractures dfnGen_output_report/networks/all_fracture_centers.pdf and one per family dfnGen_output_report/family_{family_id}/family_{family_id}_fracture_centers.pdf.
36  """
37 
38  # This keep track of the max/min radius in the entire network
39  min_radius = None
40  max_radius = None
41 
42  print("--> Plotting Fracture Radii Distributions")
43 
44  for fam in families:
45  if fam["Distribution"] != "Constant" and fam["Global Family"] > 0:
46  family_id = fam["Global Family"]
47  dist = fam["Distribution"]
48  if params["verbose"]:
49  print(
50  f"--> Working on family {family_id} which has a {dist} distribution of radii."
51  )
52 
53  # Get analytic distributions
54  if fam["Distribution"] == 'Truncated Power-Law':
55  x, pdf, cdf = tpl(fam["Alpha"], fam["Minimum Radius (m)"],
56  fam["Maximum Radius (m)"])
57 
58  elif fam["Distribution"] == 'Lognormal':
59  x, pdf, cdf = lognormal(fam["Mean"], fam["Standard Deviation"],
60  fam["Minimum Radius (m)"],
61  fam["Maximum Radius (m)"])
62 
63  elif fam["Distribution"] == 'Exponential':
64  x, pdf, cdf = exponential(fam["Lambda"],
65  fam["Minimum Radius (m)"],
66  fam["Maximum Radius (m)"])
67 
68  # Gather family radii
69  radii_all = []
70  for i in fam["fracture list - all"]:
71  radii_all.append(fractures[i]["x-radius"])
72  radii_accepted = []
73  for i in fam["fracture list - final"]:
74  radii_accepted.append(fractures[i]["x-radius"])
75 
76  #print(f"Mean: {np.mean(radii_all)}, Variance: {np.var(radii_all)}")
77  min_val = min(min(radii_all), min(radii_accepted))
78  max_val = max(max(radii_all), max(radii_accepted))
79 
80  # Get global min/max radius
81  if min_radius is None:
82  min_radius = min_val
83  else:
84  min_radius = min(min_radius, min_val)
85 
86  if max_radius is None:
87  max_radius = max_val
88  else:
89  max_radius = max(max_radius, max_val)
90 
91  fig, axs = plt.subplots(nrows=2, ncols=2, figsize=(10, 7))
92  # Create and plot histograms
93  axs[0, 0].hist(radii_all,
94  bins=num_bins,
95  color=params["all_color"],
96  range=(min_val, max_val),
97  alpha=0.7,
98  edgecolor="k",
99  linewidth=0.5,
100  density=True,
101  label=f"All Accepted Fractures")
102 
103  axs[0, 0].plot(x,
104  pdf,
105  params["analytic_color"],
106  label="Analytical PDF")
107  #axs[0, 0].legend(loc="upper right", fontsize=14)
108  axs[0, 0].grid(True, alpha=0.5)
109  axs[0, 0].set_xlabel("Fracture Radius [m]", fontsize=18)
110  axs[0, 0].set_ylabel("Density", fontsize=18)
111 
112  # Set Tick Labels. It's more complicated due to precision errors in python.
113  axs[0, 0].set_xticklabels(axs[0, 0].get_xticks().astype(int),
114  fontsize=14)
115 
116  ticks = axs[0, 0].get_yticks()
117  labels = [f"{val:0.2f}" for val in ticks]
118  axs[0, 0].set_yticklabels(labels, fontsize=14)
119 
120  # Histogram of final values
121  axs[0, 1].hist(radii_accepted,
122  bins=num_bins,
123  color=params["final_color"],
124  range=(min_val, max_val),
125  alpha=0.7,
126  edgecolor="k",
127  linewidth=0.5,
128  density=True,
129  label=f"Fractures in the Connected Network")
130 
131  axs[0, 1].plot(x,
132  pdf,
133  params["analytic_color"],
134  label="Analytical PDF")
135  #axs[0, 1].legend(loc="upper right", fontsize=14)
136  axs[0, 1].grid(True, alpha=0.5)
137  axs[0, 1].set_xlabel("Fracture Radius [m]", fontsize=18)
138  # Set Tick Labels. It's more complicated due to precision errors in python.
139  axs[0, 1].set_xticklabels(axs[0, 1].get_xticks().astype(int),
140  fontsize=14)
141 
142  ticks = axs[0, 1].get_yticks()
143  labels = [f"{val:0.2f}" for val in ticks]
144  axs[0, 1].set_yticklabels(labels, fontsize=14)
145 
146  # ECDF
147  axs[1, 0].plot(x,
148  cdf,
149  params["analytic_color"],
150  label="Analytical CDF")
151  y, ecdf = create_ecdf(radii_all)
152  axs[1, 0].plot(y,
153  ecdf,
154  label="All Fractures CDF",
155  color=params["all_color"])
156  y, ecdf = create_ecdf(radii_accepted)
157  axs[1, 0].plot(y,
158  ecdf,
159  label="Final Fractures CDF",
160  color=params["final_color"])
161 
162  #axs[1, 0].legend(loc="lower right", fontsize=14)
163  axs[1, 0].grid(True, alpha=0.5)
164  axs[1, 0].set_xlabel("Fracture Radius [m]", fontsize=18)
165  axs[1, 0].set_ylabel("Cumulative Density", fontsize=18)
166  # Set Tick Labels
167  axs[1, 0].set_xticklabels(axs[1, 0].get_xticks().astype(int),
168  fontsize=14)
169  ticks = axs[1, 0].get_yticks()
170  labels = [f"{val:0.2f}" for val in ticks]
171  axs[1, 0].set_yticklabels(labels, fontsize=14)
172 
173 
174  axs[1, 1].set_frame_on(False)
175  axs[1, 1].get_xaxis().set_visible(False)
176  axs[1, 1].get_yaxis().set_visible(False)
177  axs[1, 1].plot(0,
178  0,
179  params["all_color"],
180  linewidth=4,
181  label="All Accepted Fractures")
182  axs[1, 1].plot(0,
183  0,
184  params["final_color"],
185  linewidth=4,
186  label="Fractures in the Connected Network")
187  axs[1, 1].plot(0,
188  0,
189  params["analytic_color"],
190  linewidth=4,
191  label="Analytic Value")
192  axs[1, 1].legend(loc="center", fontsize=18, frameon=False)
193 
194  fig.tight_layout()
195  if family_id > 0:
196  fileout = f"family_{family_id}_radii.png"
197  else:
198  tmp = fam["Distribution"]
199  fileout = f"family_{tmp}_radii.png"
200 
201  if params["verbose"]:
202  print(f"--> Saving File {fileout}")
203  plt.savefig(f"{params['output_dir']}/family_{family_id}/{fileout}")
204  plt.clf()
205  plt.close()
206 
207  # Plotting all fractures
208  fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(16, 8))
209  for fam in families:
210  if fam["Global Family"] > 0 and fam["Distribution"] != "Constant":
211  family_id = fam["Global Family"]
212 
213  # Gather family radii
214  radii_all = []
215  for i in fam["fracture list - all"]:
216  radii_all.append(fractures[i]["x-radius"])
217  radii_accepted = []
218  for i in fam["fracture list - final"]:
219  radii_accepted.append(fractures[i]["x-radius"])
220 
221  min_val = min(min(radii_all), min(radii_accepted))
222  max_val = max(max(radii_all), max(radii_accepted))
223 
224  # Create and plot histograms
225  axs[0].hist(radii_all,
226  bins=num_bins,
227  color=fam["color"],
228  range=(min_radius, max_radius),
229  alpha=0.7,
230  edgecolor="k",
231  linewidth=0.5,
232  label=f"Family \# {family_id}")
233 
234  # Histogram of final values
235  axs[1].hist(radii_accepted,
236  bins=num_bins,
237  color=fam["color"],
238  range=(min_radius, max_radius),
239  alpha=0.7,
240  edgecolor="k",
241  linewidth=0.5,
242  label=f"Family \# {family_id}")
243 
244  axs[0].grid(True, alpha=0.5)
245  axs[0].set_title("All Accepted Fractures", fontsize=24)
246  axs[0].set_xlabel("Fracture Radius [m]", fontsize=24)
247  axs[0].set_ylabel("Fracture Count", fontsize=24)
248  axs[0].set_xticklabels(axs[0].get_xticks().astype(int), fontsize=14)
249  axs[0].set_yticklabels(axs[0].get_yticks().astype(int), fontsize=14)
250 
251  axs[1].grid(True, alpha=0.5)
252  axs[1].set_title("Fractures in the Connected Network", fontsize=24)
253  axs[1].set_xlabel("Fracture Radius [m]", fontsize=24)
254  axs[1].set_xticklabels(axs[1].get_xticks().astype(int), fontsize=14)
255  axs[1].set_yticklabels(axs[1].get_yticks().astype(int), fontsize=14)
256  axs[1].axis([
257  axs[0].get_xlim()[0], axs[0].get_xlim()[1], axs[0].get_ylim()[0],
258  axs[0].get_ylim()[1]
259  ])
260  plt.savefig(f"{params['output_dir']}/network/network_all_radii.png")
261  plt.clf()
262  plt.close()
def plot_fracture_radii(params, families, fractures, num_bins=20)