pydfnWorks
python wrapper for dfnWorks
wells.py
Go to the documentation of this file.
1 import os
2 import numpy as np
3 import shutil
4 
5 from pydfnworks import *
6 from pydfnworks.dfnGen.meshing import mesh_dfn_helper as mh
7 
8 
9 def tag_well_in_mesh(self, wells):
10  """ Identifies nodes in a DFN for nodes the intersect a well with radius r [m]\n
11  1. Well coordinates in well["filename"] are converted to a polyline that are written into "well_{well['name']}_line.inp"\n
12  2. Well is expanded to a volume with radius well["r"] and written into the avs file well_{well["name"]}_volume.inp\n
13  3. Nodes in the DFN that intersect with the well are written into the zone file well_{well["name"]}.zone\n
14  4. If using PFLOTRAN, then an ex file is created from the well zone file\n
15 
16  Parameters
17  -----------
18  self : object
19  DFN Class
20  well: Dictionary
21  Dictionary of information about the well that contains the following attributes
22 
23  well["name"] : string
24  name of the well
25 
26  well["filename"] : string
27  filename of the well coordinates with the following format
28  x0 y0 z0\n
29  x1 y1 z1\n
30  ...\n
31  xn yn zn\n
32 
33  well["r"] : float
34  radius of the well
35  Returns
36  --------
37  None
38  Notes
39  --------
40  Wells can be a list of well dictionaries
41 
42  """
43 
44  if type(wells) is dict:
45  well = wells
46  print(
47  f"\n\n--> Identifying nodes in the DFN intersecting with a vertical well named {well['name']}."
48  )
49 
50  # 1) convert well into polyline AVS if it doesn't exist
51  if not os.path.isfile(f"well_{well['name']}_line.inp"):
52  convert_well_to_polyline_avs(well, self.h)
53 
54  # 2) expand the polyline of the well into a volume with radius r
55  expand_well(well)
56 
57  # 3) find the nodes in the well that corresponds / intersect the well
58  get_well_zone(well, self.inp_file)
59 
60  # 4) convert the zone file to ex files for PFLTORAN
61  if self.flow_solver == "PFLOTRAN":
62  self.zone2ex(zone_file=f"well_{well['name']}.zone", face='well')
63 
64  if self.flow_solver == "FEHM":
65  print(f"--> Well nodes are in well_{well['name']}.zone")
66 
67  print(f"--> Well creation for {well['name']} complete\n\n")
68 
69  if type(wells) is list:
70  for well in wells:
71  print(
72  f"\n\n--> Identifying nodes in the DFN intersecting with a vertical well named {well['name']}."
73  )
74 
75  # 1) convert well into polyline AVS if it doesn't exist
76  if not os.path.isfile(f"well_{well['name']}_line.inp"):
77  convert_well_to_polyline_avs(well, self.h)
78 
79  # 2) expand the polyline of the well into a volume with radius r
80  expand_well(well)
81 
82  # 3) find the nodes in the well that corresponds / intersect the well
83  get_well_zone(well, self.inp_file)
84 
85  # 4) convert the zone file to ex files for PFLTORAN
86  if self.flow_solver == "PFLOTRAN":
87  self.zone2ex(zone_file=f"well_{well['name']}.zone",
88  face='well')
89 
90  if self.flow_solver == "FEHM":
91  print(f"--> Well nodes are in well_{well['name']}.zone")
92 
93  print(f"--> Well creation for {well['name']} complete\n\n")
94 
95 
96 def convert_well_to_polyline_avs(well, h=None):
97  """ Identifies converts well coordinates into a polyline avs file. Distance between
98  point on the polyline are h/2 apart. Polyline is written into "well_{well['name']}_line.inp"
99 
100  Parameters
101  -----------
102  well: dictionary of information about the well. Contains the following:
103 
104  well["name"] : string
105  name of the well
106 
107  well["filename"] : string
108  filename of the well coordinates. "well_coords.dat" for example.
109  Format is :
110  x0 y0 z0
111  x1 y1 z1
112  ...
113  xn yn zn
114 
115  well["r"] : float
116  radius of the well
117  h : float
118  h parameter for meshing.
119 
120  Returns
121  --------
122  None
123 
124  Notes
125  --------
126  If flow solver is set to PFLOTRAN, the zone file dumped by LaGriT will be converted to
127  an ex file.
128 
129  """
130  print("--> Interpolating well coordinates into a polyline")
131  # If h is not provided, get it from params.txt
132  if h == None:
133  from pydfnworks.dfnGen.meshing.mesh_dfn_helper import parse_params_file
134  _, h, _, _, _ = parse_params_file(quiet=True)
135 
136  # read in well coordinates
137  pts = np.genfromtxt(f"{well['filename']}")
138  n, _ = np.shape(pts)
139 
140  # Linear interpolation of well into a polyline
141  new_pts = []
142  new_pts.append(pts[0])
143  new_idx = 0
144 
145  for i in range(1, n):
146  distance = np.linalg.norm(pts[i, :] - pts[i - 1, :])
147  if distance < h:
148  new_pts.append(pts[i, :])
149  new_idx += 1
150  else:
151  # discretized to less than h
152  m = int(np.ceil(distance / h))
153  dx = (pts[i, 0] - pts[i - 1, 0]) / m
154  dy = (pts[i, 1] - pts[i - 1, 1]) / m
155  dz = (pts[i, 2] - pts[i - 1, 2]) / m
156  for j in range(m):
157  interp = np.zeros(3)
158  interp[0] = new_pts[new_idx][0] + dx
159  interp[1] = new_pts[new_idx][1] + dy
160  interp[2] = new_pts[new_idx][2] + dz
161  new_pts.append(interp)
162  del interp
163  new_idx += 1
164 
165  print("--> Interpolating well coordinates into a polyline: Complete")
166  # Write interpolated polyline into an AVS file
167  avs_filename = f"well_{well['name']}_line.inp"
168  print(f"--> Writing polyline into avs file : {avs_filename}")
169 
170  num_pts = new_idx + 1
171  pt_digits = len(str(num_pts))
172 
173  num_elem = new_idx
174  elem_digits = len(str(num_elem))
175 
176  favs = open(avs_filename, "w")
177  favs.write(f"{num_pts}\t{num_elem}\t0\t0\t0\n")
178  for i in range(num_pts):
179  favs.write(
180  f"{i+1:0{pt_digits}d} {new_pts[i][0]:12e} {new_pts[i][1]:12e} {new_pts[i][2]:12e}\n"
181  )
182  for i in range(num_elem):
183  favs.write(f"{i+1} 1 line {i+1} {i+2}\n")
184  favs.close()
185  print(f"--> Writing polyline into avs file : {avs_filename} : Complete")
186 
187 
188 def expand_well(well):
189  """ Expands the polyline defining the well into a volume with radius r [m].
190  A sphere of points around each point is created and then connected.
191  Volume is written into the avs file well_{well["name"]}_volume.inp
192 
193  Parameters
194  -----------
195  well:
196  dictionary of information about the well. Contains the following:
197 
198  well["name"] : string
199  name of the well
200 
201  well["filename"] : string
202  filename of the well coordinates. "well_coords.dat" for example.
203  Format is :
204  x0 y0 z0
205  x1 y1 z1
206  ...
207  xn yn zn
208 
209  well["r"] : float
210  radius of the well
211 
212  Returns
213  --------
214  None
215 
216  Notes
217  --------
218  Mesh of the well is written into the avs file {well["name"][:-4]}.inp
219 
220  """
221 
222  print("--> Expanding well into a volume.")
223 
224  r = well["r"]
225 
226  convert_well_to_polyline_avs(well, r)
227 
228  angle_r = r / np.sqrt(2)
229 
230  lagrit_script = f"""
231  ## read in polyline of well
232  read / well_{well['name']}_line.inp / mo_line
233  cmo / printatt / mo_line / -xyz- / minmax
234 
235  ## expand every point in the polyline into a discrete sphere
236  ## with radius r
237  cmo / create / mo_well / / / tet
238  copypts / mo_well / mo_line
239 
240  cmo / create / mo_tmp / / / line
241  copypts / mo_tmp / mo_line
242  cmo / select / mo_tmp
243  trans / 1 0 0 / 0. 0. 0. / {well["r"]} 0 0
244  copypts / mo_well / mo_tmp
245  cmo / delete / mo_tmp
246 
247  cmo / create / mo_tmp / / / line
248  copypts / mo_tmp / mo_line
249  cmo / select / mo_tmp
250  trans / 1 0 0 / 0. 0. 0. / -{well["r"]} 0 0
251  copypts / mo_well / mo_tmp
252  cmo / delete / mo_tmp
253 
254  cmo / create / mo_tmp / / / line
255  copypts / mo_tmp / mo_line
256  cmo / select / mo_tmp
257  trans / 1 0 0 / 0. 0. 0. / 0 {well["r"]} 0
258  copypts / mo_well / mo_tmp
259  cmo / delete / mo_tmp
260 
261  cmo / create / mo_tmp / / / line
262  copypts / mo_tmp / mo_line
263  cmo / select / mo_tmp
264  trans / 1 0 0 / 0. 0. 0. / 0 -{well["r"]} 0
265  copypts / mo_well / mo_tmp
266  cmo / delete / mo_tmp
267 
268  cmo / create / mo_tmp / / / line
269  copypts / mo_tmp / mo_line
270  cmo / select / mo_tmp
271  trans / 1 0 0 / 0. 0. 0. / 0 0 {well["r"]}
272  copypts / mo_well / mo_tmp
273  cmo / delete / mo_tmp
274 
275  cmo / create / mo_tmp / / / line
276  copypts / mo_tmp / mo_line
277  cmo / select / mo_tmp
278  trans / 1 0 0 / 0. 0. 0. / 0 0 -{well["r"]}
279  copypts / mo_well / mo_tmp
280  cmo / delete / mo_tmp
281 
282  cmo / create / mo_tmp / / / line
283  copypts / mo_tmp / mo_line
284  cmo / select / mo_tmp
285  trans / 1 0 0 / 0. 0. 0. / {angle_r} {angle_r} 0
286  copypts / mo_well / mo_tmp
287  cmo / delete / mo_tmp
288 
289  cmo / create / mo_tmp / / / line
290  copypts / mo_tmp / mo_line
291  cmo / select / mo_tmp
292  trans / 1 0 0 / 0. 0. 0. / {angle_r} -{angle_r} 0
293  copypts / mo_well / mo_tmp
294  cmo / delete / mo_tmp
295 
296  cmo / create / mo_tmp / / / line
297  copypts / mo_tmp / mo_line
298  cmo / select / mo_tmp
299  trans / 1 0 0 / 0. 0. 0. / -{angle_r} {angle_r} 0
300  copypts / mo_well / mo_tmp
301  cmo / delete / mo_tmp
302 
303  cmo / create / mo_tmp / / / line
304  copypts / mo_tmp / mo_line
305  cmo / select / mo_tmp
306  trans / 1 0 0 / 0. 0. 0. / -{angle_r} -{angle_r} 0
307  copypts / mo_well / mo_tmp
308  cmo / delete / mo_tmp
309 
310  cmo / create / mo_tmp / / / line
311  copypts / mo_tmp / mo_line
312  cmo / select / mo_tmp
313  trans / 1 0 0 / 0. 0. 0. / 0 {angle_r} {angle_r}
314  copypts / mo_well / mo_tmp
315  cmo / delete / mo_tmp
316 
317  cmo / create / mo_tmp / / / line
318  copypts / mo_tmp / mo_line
319  cmo / select / mo_tmp
320  trans / 1 0 0 / 0. 0. 0. / 0 {angle_r} -{angle_r}
321  copypts / mo_well / mo_tmp
322  cmo / delete / mo_tmp
323 
324  cmo / create / mo_tmp / / / line
325  copypts / mo_tmp / mo_line
326  cmo / select / mo_tmp
327  trans / 1 0 0 / 0. 0. 0. / 0 -{angle_r} {angle_r}
328  copypts / mo_well / mo_tmp
329  cmo / delete / mo_tmp
330 
331  cmo / create / mo_tmp / / / line
332  copypts / mo_tmp / mo_line
333  cmo / select / mo_tmp
334  trans / 1 0 0 / 0. 0. 0. / 0 -{angle_r} -{angle_r}
335  copypts / mo_well / mo_tmp
336  cmo / delete / mo_tmp
337 
338  cmo / create / mo_tmp / / / line
339  copypts / mo_tmp / mo_line
340  cmo / select / mo_tmp
341  trans / 1 0 0 / 0. 0. 0. / {angle_r} 0 {angle_r}
342  copypts / mo_well / mo_tmp
343  cmo / delete / mo_tmp
344 
345  cmo / create / mo_tmp / / / line
346  copypts / mo_tmp / mo_line
347  cmo / select / mo_tmp
348  trans / 1 0 0 / 0. 0. 0. / {angle_r} 0 -{angle_r}
349  copypts / mo_well / mo_tmp
350  cmo / delete / mo_tmp
351 
352  cmo / create / mo_tmp / / / line
353  copypts / mo_tmp / mo_line
354  cmo / select / mo_tmp
355  trans / 1 0 0 / 0. 0. 0. / -{angle_r} 0 {angle_r}
356  copypts / mo_well / mo_tmp
357  cmo / delete / mo_tmp
358 
359  cmo / create / mo_tmp / / / line
360  copypts / mo_tmp / mo_line
361  cmo / select / mo_tmp
362  trans / 1 0 0 / 0. 0. 0. / -{angle_r} 0 -{angle_r}
363  copypts / mo_well / mo_tmp
364  cmo / delete / mo_tmp
365 
366  ## Could add more permutations, but this looks good enough right now
367  ## JDH 9 Sept. 2020
368 
369  ################## DEBUG ###########################
370  # dump / well_pts.inp / mo_well
371  ################## DEBUG ###########################
372 
373  # filter out point that are too close
374  cmo / select / mo_well
375  filter / 1 0 0 / {0.1*r}
376  rmpoint / compress
377  cmo / setatt / mo_well / imt / 1 0 0 / 1
378 
379  # connect the point cloud and make a volume mesh
380  connect / delaunay
381  resetpts / itp
382 
383  ################## DEBUG ###########################
384  # dump / well_{well["name"]}.inp / mo_well
385  ################## DEBUG ###########################
386 
387  # add edge_max attribute and remove elements with big edges
388  quality / edge_max / y
389  eltset /big/ edgemax / gt / {np.sqrt(2)*r}
390  cmo / setatt / mo_well / itetclr / eltset, get, big / 2
391  rmmat / 2
392  rmpoint / compress
393 
394  dump / well_{well["name"]}_volume.inp / mo_well
395 
396  ################## DEBUG ###########################
397  # # extract the surface of the well
398  # # This is done to remove internal points and reduce
399  # # the total number of elements in the mesh
400  # # This speeds up the intersection checking later on
401  # # I couldn't get this to work in a robust way.
402  # # There were some weird LaGriT errors if I deleted
403  # # mesh object.
404  # # Works if we stop before this, but I'm leaving it to
405  # # revisit if need be.
406  # # JDH 10/9/2020
407 
408  # extract / surfmesh / 1,0,0 /mo_shell / mo_well
409  # cmo / select / mo_shell
410  # cmo / delete / mo_well
411 
412  # ################## DEBUG ###########################
413  # dump / well_{well["name"]}_shell.inp / mo_shell
414  # ################## DEBUG ###########################
415 
416  # # Copy the surface of the well into a tet mesh
417  # cmo / create / mo_well2 / / / tet
418  # copypts / mo_well2 / mo_shell
419  # cmo / select / mo_well2
420  # # cmo / delete / mo_shell
421 
422  # # filter out point that are too close
423  # filter / 1 0 0 / {0.1*r}
424  # rmpoint / compress
425  # cmo / setatt / mo_well2 / imt / 1 0 0 / 1
426 
427  # # connect the point cloud and make a volume mesh
428  # connect / delaunay
429  # resetpts / itp
430 
431  # # add edge_max attribute and remove elements with big edges
432  # quality / edge_max / y
433  # eltset /big/ edgemax / gt / {np.sqrt(2)*r}
434  # cmo / setatt / mo_well2 / itetclr / eltset, get, big / 2
435  # rmmat / 2
436  # rmpoint / compress
437 
438  # # write out final well mesh
439  # dump / well_{well["name"]}_volume.inp / mo_well2
440 
441  finish
442 
443  """
444  # Write LaGriT commands to file
445  with open(f"expand_well_{well['name']}.lgi", "w") as fp:
446  fp.write(lagrit_script)
447  fp.flush()
448 
449  # Execute LaGriT
450  mh.run_lagrit_script(f"expand_well_{well['name']}.lgi",
451  output_file=f"expand_well_{well['name']}.out",
452  quiet=False)
453 
454  print("--> Expanding well complete.")
455 
456 
457 def get_well_zone(well, inp_file):
458  """Identifies nodes in a DFN for nodes the intersect a well with radius r [m]
459  First, all elements that intersect the well are identified.
460  Second, all nodes of those elements are tagged.
461  Third, that collection of nodes are dumped as a zone file (well_{well["name"]}.zone)
462 
463  Parameters
464  -----------
465  self : object
466  DFN Class
467  well:
468  dictionary of information about the well. Contains the following:
469 
470  well["name"] : string
471  name of the well
472 
473  well["filename"] : string
474  filename of the well coordinates. "well_coords.dat" for example.
475  File format:
476  x0 y0 z0
477  x1 y1 z1
478  ...
479  xn yn zn
480 
481  well["r"] : float
482  radius of the well
483 
484  Returns
485  --------
486  None
487 
488  Notes
489  --------
490  None
491 
492  """
493  # # if the well has not been converted to AVS, do that first
494  # if not os.path.isfile(f"well_{well['name']}_line.inp"):
495  # convert_well_to_polyline_avs(well,h)
496  # # if the well has not been expanded
497  # if not os.path.isfile(f"well_{well['name']}_volume.inp"):
498  # expand_well(well)
499 
500  lagrit_script = f"""
501 
502 # read in well volume
503 read / well_{well["name"]}_volume.inp / mo_well
504 
505 # read in DFN
506 read / {inp_file} / mo_dfn
507 
508 # find intersecting cells
509 cmo / select / mo_dfn
510 intersect_elements / mo_dfn / mo_well / well_{well["name"]}_inter
511 
512 eltset / ewell / well_{well["name"]}_inter / gt / 0
513 
514 # dump dfn mesh with intersections tagged
515 #dump / avs / {inp_file[:-3]}_tagged.inp / mo_dfn
516 
517 # gather nodes of intersecting cells
518 pset / well_{well["name"]} / eltset / ewell
519 
520 # dump nodes from intersecting cells
521 pset / well_{well["name"]} / zone / well_{well["name"]}.zone
522 
523 finish
524 
525 """
526  # Write LaGriT commands to file
527  with open(f"get_well_{well['name']}_zone.lgi", "w") as fp:
528  fp.write(lagrit_script)
529  fp.flush()
530 
531  # Execute LaGriT
532  mh.run_lagrit_script(f"get_well_{well['name']}_zone.lgi",
533  output_file=f"create_well_{well['name']}.out",
534  quiet=False)
535 
536  with open(f"well_{well['name']}.zone", "r") as fp:
537  number_of_nodes = int(fp.readlines()[3])
538  if number_of_nodes > 0:
539  print(f"--> There are {number_of_nodes} nodes in the well zone")
540  else:
541  print("--> WARNING!!! The well did not intersect the DFN!!!")
542 
543 
544 def find_well_intersection_points(self, wells):
545  """ Identifies points on a DFN where the well intersects the network.
546  These points are used in meshing the network to have higher resolution in the mesh in these points.
547  Calls a sub-routine run_find_well_intersection_points.
548 
549  Parameters
550  -----------
551  self : object
552  DFN Class
553  well:
554  dictionary of information about the well. Contains the following:
555 
556  well["name"] : string
557  name of the well
558 
559  well["filename"] : string
560  filename of the well coordinates. "well_coords.dat" for example.
561  Format is :
562  x0 y0 z0
563  x1 y1 z1
564  ...
565  xn yn zn
566 
567  well["r"] : float
568  radius of the well
569 
570  Returns
571  --------
572  None
573 
574  Notes
575  --------
576  Wells can be a list of well dictionaries.
577  Calls the subroutine run_find_well_intersection_points to remove redundant code.
578 
579  """
580  # check for reduced mesh, if it doesn't exists, make it
581  print("--> Checking for reduced_mesh.inp")
582  if not os.path.isfile("reduced_mesh.inp"):
583  print("--> reduced_mesh.inp not found. Creating it now.")
584  self.mesh_network(visual_mode=True)
585  else:
586  print("--> reduced_mesh.inp found. Moving on.")
587 
588  # if using a single well
589  if type(wells) is dict:
590  run_find_well_intersection_points(wells, self.h)
591  # using a list of wells, loop over them.
592  elif type(wells) is list:
593  for well in wells:
594  run_find_well_intersection_points(well, self.h)
595 
596  # Run cross check
597  cross_check_pts(self.h)
598 
599 
600 def run_find_well_intersection_points(well, h):
601  """ Runs the workflow for finding the point of intersection of the DFN with the well.
602 
603  Parameters
604  -----------
605  well:
606  dictionary of information about the well. Contains the following:
607 
608  well["name"] : string
609  name of the well
610 
611  well["filename"] : string
612  filename of the well coordinates. "well_coords.dat" for example.
613  Format is :
614  x0 y0 z0
615  x1 y1 z1
616  ...
617  xn yn zn
618 
619  well["r"] : float
620  radius of the well
621  h : float
622  Minimum h length scale in the network
623 
624  Returns
625  --------
626  None
627 
628  Notes
629  --------
630  This function was designed to minimize redundancy in the workflow. It's called
631  within find_well_intersection_points()
632  """
633 
634  print(f"\n\n--> Working on well {well['name']}")
635 
636  # 1) convert well into polyline AVS
637  if not os.path.isfile(f"well_{well['name']}_line.inp"):
638  convert_well_to_polyline_avs(well, h)
639 
640  # run LaGriT scripts to dump information
641  find_segments(well)
642 
643  well_point_of_intersection(well)
644 
645 
646 def find_segments(well):
647  """ LaGriT script to identify the points of intersection between the DFN and the well.
648 
649  Parameters
650  -----------
651  well:
652  dictionary of information about the well. Contains the following:
653 
654  well["name"] : string
655  name of the well
656 
657  well["filename"] : string
658  filename of the well coordinates. "well_coords.dat" for example.
659  Format is :
660  x0 y0 z0
661  x1 y1 z1
662  ...
663  xn yn zn
664 
665  well["r"] : float
666  radius of the well
667 
668  Returns
669  --------
670  None
671 
672  Notes
673  --------
674  points of intersection are written into the avs file well_{well['name']}_intersect.inp
675 
676  OUTPUT: well_segments_intersect.inp is a subset of the well line segments that intersect fractures.
677  The segments are tagged so itetclr and imt are set to the value of the fracture they intersect.
678  """
679 
680  lagrit_script = f"""
681 #
682 # OUTPUT: intersected_fracture.list tells you the list of fractures intersected by the well.
683 # OUTPUT: well_segments_intersect.inp is a subset of the well line segments that intersect fractures.
684 # The segments are tagged so itetclr and imt are set to the value of the fracture they intersect.
685 #
686 define / INPUT_DFN / reduced_mesh.inp
687 define / INPUT_WELL / well_{well['name']}_line.inp
688 define / OUTPUT_WELL_SEGMENTS / well_{well['name']}_intersect.inp
689 # define / OUTPUT_FRACTURE_LIST / {well['name']}_fracture.list
690 #
691 read / avs / INPUT_DFN / mo_tri
692 read / avs / INPUT_WELL / mo_line
693 #
694 # Find the triangles of the DFN mesh that intersect the well lines.
695 # Get rid of all the non-intersecting triangles.
696 #
697 intersect_elements / mo_tri / mo_line / if_intersect
698 eltset / e_not_intersect / if_intersect / eq / 0
699 rmpoint / element / eltset get e_not_intersect
700 rmpoint / compress
701 cmo / DELATT / mo_tri / if_intersect
702 #
703 # dump / avs / reduced_reduced_mesh.inp / mo_tri
704 cmo / addatt / mo_tri / id_fracture / vint / scalar / nelements
705 cmo / copyatt / mo_tri / mo_tri / id_fracture / itetclr
706 # dump / avs / OUTPUT_FRACTURE_LIST / mo_tri / 0 0 0 1
707 #
708 # Find the segments of the well (line object) that intersect the fracture planes (triangles)
709 #
710 intersect_elements / mo_line / mo_tri / if_intersect
711 eltset / e_not_intersect / if_intersect / eq / 0
712 rmpoint / element / eltset get e_not_intersect
713 rmpoint / compress
714 cmo / DELATT / mo_line / if_intersect
715 # BEGIN DEBUG
716 # dump / avs / OUTPUT_WELL_SEGMENTS / mo_line
717 # END DEBUG
718 #
719 # Reduce the size of the triangles so interpolation works.
720 #
721 cmo / select / mo_tri
722 # Refine 2**7 128
723 refine2d
724 intersect_elements / mo_tri / mo_line / if_intersect
725 eltset / e_not_intersect / if_intersect / eq / 0
726 rmpoint / element / eltset get e_not_intersect
727 rmpoint / compress
728 cmo / DELATT / mo_tri / if_intersect
729 
730 refine2d
731 intersect_elements / mo_tri / mo_line / if_intersect
732 eltset / e_not_intersect / if_intersect / eq / 0
733 rmpoint / element / eltset get e_not_intersect
734 rmpoint / compress
735 cmo / DELATT / mo_tri / if_intersect
736 
737 refine2d
738 intersect_elements / mo_tri / mo_line / if_intersect
739 eltset / e_not_intersect / if_intersect / eq / 0
740 rmpoint / element / eltset get e_not_intersect
741 rmpoint / compress
742 cmo / DELATT / mo_tri / if_intersect
743 
744 refine2d
745 intersect_elements / mo_tri / mo_line / if_intersect
746 eltset / e_not_intersect / if_intersect / eq / 0
747 rmpoint / element / eltset get e_not_intersect
748 rmpoint / compress
749 cmo / DELATT / mo_tri / if_intersect
750 
751 refine2d
752 intersect_elements / mo_tri / mo_line / if_intersect
753 eltset / e_not_intersect / if_intersect / eq / 0
754 rmpoint / element / eltset get e_not_intersect
755 rmpoint / compress
756 cmo / DELATT / mo_tri / if_intersect
757 
758 refine2d
759 intersect_elements / mo_tri / mo_line / if_intersect
760 eltset / e_not_intersect / if_intersect / eq / 0
761 rmpoint / element / eltset get e_not_intersect
762 rmpoint / compress
763 cmo / DELATT / mo_tri / if_intersect
764 
765 refine2d
766 intersect_elements / mo_tri / mo_line / if_intersect
767 eltset / e_not_intersect / if_intersect / eq / 0
768 rmpoint / element / eltset get e_not_intersect
769 rmpoint / compress
770 cmo / DELATT / mo_tri / if_intersect
771 
772 # BEGIN DEBUG
773 # dump / avs / tmp_refine.inp / mo_tri
774 # END DEBUG
775 
776 interpolate / voronoi / mo_line itetclr / 1 0 0 / mo_tri imt
777 interpolate / voronoi / mo_line imt / 1 0 0 / mo_tri imt
778 
779 cmo / modatt / mo_line / itp / ioflag / l
780 cmo / modatt / mo_line / icr / ioflag / l
781 cmo / modatt / mo_line / isn / ioflag / l
782 
783 dump / avs / OUTPUT_WELL_SEGMENTS / mo_line
784 
785 finish
786 
787 
788 """
789  # Write LaGriT commands to file
790  with open(f"find_well_{well['name']}_segment.lgi", "w") as fp:
791  fp.write(lagrit_script)
792  fp.flush()
793 
794  # Execute LaGriT
795  mh.run_lagrit_script(f"find_well_{well['name']}_segment.lgi",
796  output_file=f"find_well_{well['name']}_segment.out",
797  quiet=False)
798 
799 
800 def well_point_of_intersection(well):
801  """ Takes the well points found using find_segments and projects the points onto the fracture plane. These points are written into well_points.dat file. During meshing, these points are read in and a higher resolution mesh is created near by them. well_points.dat has the format
802 
803  fracture_id x y z
804  ...
805 
806  for every intersection point.
807 
808  Parameters
809  -----------
810  well:
811  dictionary of information about the well. Contains the following:
812 
813  well["name"] : string
814  name of the well
815 
816  well["filename"] : string
817  filename of the well coordinates. "well_coords.dat" for example.
818  Format is :
819  x0 y0 z0
820  x1 y1 z1
821  ...
822  xn yn zn
823 
824  well["r"] : float
825  radius of the well
826 
827  Returns
828  --------
829  None
830 
831  Notes
832  --------
833 
834  """
835 
836  print(f"--> Finding well points on DFN for {well['name']}")
837  # create file to keep well points if it doesn't exist. Otherwise set to append.
838  if not os.path.isfile("well_points.dat"):
839  fwell = open("well_points.dat", "w")
840  fwell.write("fracture_id x y z\n")
841  else:
842  fwell = open("well_points.dat", "a")
843 
844  well_line_file = f"well_{well['name']}_intersect.inp"
845 
846  pts, elems, fracture_list = get_segments(well_line_file)
847 
848  if len(fracture_list) == 0:
849  print(
850  f"\n--> WARNING!!! The well {well['name']} did not intersect the DFN!!!\n"
851  )
852 
853  for elem in elems: # Parameterize the line center of the well
854  l0 = np.zeros(3)
855  l0[0] = pts[elem["pt1"] - 1]["x"]
856  l0[1] = pts[elem["pt1"] - 1]["y"]
857  l0[2] = pts[elem["pt1"] - 1]["z"]
858  l1 = np.zeros(3)
859  l1[0] = pts[elem["pt2"] - 1]["x"]
860  l1[1] = pts[elem["pt2"] - 1]["y"]
861  l1[2] = pts[elem["pt2"] - 1]["z"]
862  l = l1 - l0
863 
864  f = elem["frac"]
865  # get the plane on which the fracture lies
866  n = get_normal(f)
867  p0 = get_center(f)
868  R = rotation_matrix(n, [0, 0, 1])
869 
870  # find the point of intersection between the well line and the plane
871  d = np.dot((p0 - l0), n) / (np.dot(l, n))
872  p = l0 + l * d
873  v = rotate_point(p, R)
874  fwell.write(f"{f} {v[0]} {v[1]} {v[2]}\n")
875  fwell.close()
876 
877 
878 def cross_check_pts(h):
879  """ Sometimes multiple points of intersection are identified on the same fracture. This can occur if the discretized well has points close to the fracture plane. This function walks through well_points.dat and removes duplicate points that are within h of one another and on the same fracture plane.
880 
881 
882  Parameters
883  -----------
884  h : float
885  Minimum length scale in the network.
886 
887  Returns
888  --------
889  None
890 
891  Notes
892  --------
893  None
894  """
895 
896  print("\n--> Cross Checking well points")
897  pts = np.genfromtxt("well_points.dat", skip_header=1)
898  num_pts, _ = np.shape(pts)
899 
900  # Walk through well points and see if they are too close together,
901  # This can happen due to machine precision in LaGriT.
902  # We only keep 1 point per fracture per well.
903  remove_idx = []
904  for i in range(num_pts):
905  fracture_number = pts[i, 0]
906  for j in range(i + 1, num_pts):
907  # Check if points are on the same fracture to reduce number of distance checks
908  if fracture_number == pts[j, 0]:
909  dist = abs(pts[i, 1] - pts[j, 1]) + abs(
910  pts[i, 2] - pts[j, 2]) + abs(pts[i, 3] - pts[j, 3])
911  # if the points are closure that h/2, mark one to be removed.
912  if dist < h / 2:
913  remove_idx.append(j)
914 
915  # Find the id of those points to keep
916  keep_pt_indes = list(set(range(num_pts)) - set(remove_idx))
917  # keep only those points
918  pts = pts[keep_pt_indes]
919  num_pts = len(pts)
920  # write to file
921  # os.remove("well_points.dat")
922  with open("well_points.dat", "w") as fwell:
923  fwell.write("fracture_id x y z\n")
924  for i in range(num_pts):
925  fwell.write(f"{int(pts[i,0])} {pts[i,1]} {pts[i,2]} {pts[i,3]}\n")
926  print("--> Cross Checking Complete")
927 
928 
929 def get_segments(well_line_file):
930  """ Parses well_line_file (avs) to get point information, element information, and list of fractures that intersect the well.
931 
932  Parameters
933  -----------
934  well_line_file : string
935  filename of well_line_file written by find_segments()
936 
937  Returns
938  --------
939  pts : list
940  list of dictionaries of intersection points. dictionary contains fracture id and x,y,z coordinates
941  elems: list of dictionaries
942  Information about elements of the discretized well that intersect the DFN
943  fracture_list : list
944  list of fractures that the well intersects
945 
946  Notes
947  --------
948  None
949  """
950 
951  with open(well_line_file, "r") as fp:
952  header = fp.readline().split()
953  num_pts = int(header[0])
954  num_elem = int(header[1])
955  pts = []
956  for i in range(num_pts):
957  pt = fp.readline().split()
958  tmp = {"id": None, "x": None, "y": None, "z": None}
959  tmp["id"] = int(pt[0])
960  tmp["x"] = float(pt[1])
961  tmp["y"] = float(pt[2])
962  tmp["z"] = float(pt[3])
963  pts.append(tmp)
964  elems = []
965  for i in range(num_elem):
966  elem = fp.readline().split()
967  tmp = {"pt1": None, "pt2": None, "frac": None}
968  tmp["pt1"] = int(elem[3])
969  tmp["pt2"] = int(elem[4])
970  tmp["frac"] = int(elem[1])
971  elems.append(tmp)
972  # get fracture list
973  fracture_list = []
974  for i in range(num_elem):
975  if not elems[i]["frac"] in fracture_list:
976  fracture_list.append(elems[i]["frac"])
977 
978  return pts, elems, fracture_list
979 
980 
981 def get_normal(fracture_id):
982  """ Returns Normal vector of a fracture
983 
984  Parameters
985  -----------
986  fracture_id : int
987  fracture number
988 
989  Returns
990  --------
991  normal : numpy array
992  normal vector of a fracture
993 
994  Notes
995  --------
996  None
997  """
998 
999  normals = np.genfromtxt("normal_vectors.dat")
1000  return normals[fracture_id - 1, :]
1001 
1002 
1003 def get_center(fracture_id):
1004  """ Returns center of a fracture
1005 
1006  Parameters
1007  -----------
1008  fracture_id : int
1009  fracture number
1010 
1011  Returns
1012  --------
1013  points : numpy array
1014  x,y,z coordinates of a fracture
1015 
1016  Notes
1017  --------
1018  None
1019  """
1020  with open('translations.dat') as old, open('points.dat', 'w') as new:
1021  old.readline()
1022  for line in old:
1023  if not 'R' in line:
1024  new.write(line)
1025  points = np.genfromtxt('points.dat', skip_header=0, delimiter=' ')
1026  return points[fracture_id - 1, :]
1027 
1028 
1029 def rotation_matrix(normalA, normalB):
1030  """ Create a Rotation matrix to transform normal vector A to normal vector B
1031 
1032  Parameters
1033  -----------
1034  normalA : numpy array
1035  normal vector
1036  normalB : numpy array
1037  normal vector
1038 
1039 
1040  Returns
1041  --------
1042  R : numpy array
1043  Rotation matrix
1044 
1045  Notes
1046  --------
1047  None
1048  """
1049  # Check if normals are the same.
1050  comparison = normalA == normalB
1051  equal_arrays = comparison.all()
1052 
1053  # If they are equal, Return the Identity Matrix
1054  if equal_arrays:
1055  R = np.zeros(9)
1056  R[0] = 1
1057  R[1] = 0
1058  R[2] = 0
1059  R[3] = 0
1060  R[4] = 1
1061  R[5] = 0
1062  R[6] = 0
1063  R[7] = 0
1064  R[8] = 1
1065  # If they are not equal, construct and return a Rotation Matrix
1066  else:
1067 
1068  xProd = np.cross(normalA, normalB)
1069  sin = np.sqrt(xProd[0] * xProd[0] + xProd[1] * xProd[1] +
1070  xProd[2] * xProd[2])
1071  cos = np.dot(normalA, normalB)
1072  v = np.zeros(9)
1073  v = [
1074  0, -xProd[2], xProd[1], xProd[2], 0, -xProd[0], -xProd[1],
1075  xProd[0], 0
1076  ]
1077  scalar = (1.0 - cos) / (sin * sin)
1078  vSquared = np.zeros(9)
1079  vSquared[0] = (v[0] * v[0] + v[1] * v[3] + v[2] * v[6]) * scalar
1080  vSquared[1] = (v[0] * v[1] + v[1] * v[4] + v[2] * v[7]) * scalar
1081  vSquared[2] = (v[0] * v[2] + v[1] * v[5] + v[2] * v[8]) * scalar
1082  vSquared[3] = (v[3] * v[0] + v[4] * v[3] + v[5] * v[6]) * scalar
1083  vSquared[4] = (v[3] * v[1] + v[4] * v[4] + v[5] * v[7]) * scalar
1084  vSquared[5] = (v[3] * v[2] + v[4] * v[5] + v[5] * v[8]) * scalar
1085  vSquared[6] = (v[6] * v[0] + v[7] * v[3] + v[8] * v[6]) * scalar
1086  vSquared[7] = (v[6] * v[1] + v[7] * v[4] + v[8] * v[7]) * scalar
1087  vSquared[8] = (v[6] * v[2] + v[7] * v[5] + v[8] * v[8]) * scalar
1088  R = np.zeros(9)
1089  R[0] = 1 + v[0] + vSquared[0]
1090  R[1] = 0 + v[1] + vSquared[1]
1091  R[2] = 0 + v[2] + vSquared[2]
1092  R[3] = 0 + v[3] + vSquared[3]
1093  R[4] = 1 + v[4] + vSquared[4]
1094  R[5] = 0 + v[5] + vSquared[5]
1095  R[6] = 0 + v[6] + vSquared[6]
1096  R[7] = 0 + v[7] + vSquared[7]
1097  R[8] = 1 + v[8] + vSquared[8]
1098 
1099  return R
1100 
1101 
1102 def rotate_point(p, R):
1103  """ Apply Rotation matrix R to the point p
1104 
1105  Parameters
1106  -----------
1107  p : numpy array
1108  point in 3D space
1109  R : numpy array
1110  Rotation matrix
1111  Returns
1112  --------
1113  v : numpy array
1114  The point p with the rotation matrix applied
1115 
1116  Notes
1117  --------
1118  None
1119  """
1120  v = np.zeros(3)
1121  v[0] = p[0] * R[0] + p[1] * R[1] + p[2] * R[2]
1122  v[1] = p[0] * R[3] + p[1] * R[4] + p[2] * R[5]
1123  v[2] = p[0] * R[6] + p[1] * R[7] + p[2] * R[8]
1124  return v
1125 
1126 
1127 def cleanup_wells(self, wells):
1128  """ Moves working files created while making wells into well_data directory
1129 
1130  Parameters
1131  -----------
1132  self : object
1133  DFN Class
1134  well:
1135  dictionary of information about the well. Contains the following:
1136 
1137  well["name"] : string
1138  name of the well
1139 
1140  well["filename"] : string
1141  filename of the well coordinates. "well_coords.dat" for example.
1142  Format is :
1143  x0 y0 z0
1144  x1 y1 z1
1145  ...
1146  xn yn zn
1147 
1148  well["r"] : float
1149  radius of the well
1150 
1151  Returns
1152  --------
1153  None
1154 
1155  Notes
1156  --------
1157  Wells can be a list of well dictionaries
1158  """
1159 
1160  if not os.path.isdir("well_data"):
1161  os.mkdir("well_data")
1162 
1163  files = ["well_{0}_line.inp", "expand_well_{0}.lgi", \
1164  "well_{0}_volume.inp","expand_well_{0}.out",\
1165  "get_well_{0}_zone.lgi", "create_well_{0}.out",\
1166  "well_{0}_intersect.inp"]
1167 
1168  if type(wells) is dict:
1169  well = wells
1170  for file in files:
1171  shutil.move(file.format(well['name']),
1172  "well_data/" + file.format(well['name']))
1173 
1174  if type(wells) is list:
1175  for well in wells:
1176  for file in files:
1177  shutil.move(file.format(well['name']),
1178  "well_data/" + file.format(well['name']))
1179 
1180 
1181 def combine_well_boundary_zones(self, wells):
1182  """ Processes zone files for particle tracking. All zone files are combined into allboundaries.zone
1183 
1184  Parameters
1185  ----------
1186  None
1187 
1188  Returns
1189  -------
1190  None
1191 
1192  Notes
1193  -----
1194  None
1195  """
1196  # If there is only 1 well, make a symbolic link
1197  if type(wells) is dict:
1198  os.symlink(f"well_{well['name']}.zone", "well_nodes.zone")
1199 
1200  if type(wells) is list:
1201  number_of_wells = len(wells)
1202  fall = open("well_nodes.zone", "w")
1203  for index, well in enumerate(wells):
1204  if index == 0:
1205  print(f"Working on well {well['name']}")
1206  fzone = open(f"well_{well['name']}.zone", "r")
1207  lines = fzone.readlines()
1208  lines = lines[:-2]
1209  fall.writelines(lines)
1210  if index > 0 and index < number_of_wells - 1:
1211  print(f"Working on well {well['name']}")
1212  fzone = open(f"well_{well['name']}.zone", "r")
1213  lines = fzone.readlines()
1214  lines = lines[1:-2]
1215  lines[0] = f"{index+1:06d}\t\t{well['name']}\n"
1216  fzone.close()
1217  fall.writelines(lines)
1218  if index == number_of_wells - 1:
1219  print(f"Working on well {well['name']}")
1220  fzone = open(f"well_{well['name']}.zone", "r")
1221  lines = fzone.readlines()
1222  lines = lines[1:]
1223  lines[0] = f"{index+1:06d}\t\t{well['name']}\n"
1224  fzone.close()
1225  fall.writelines(lines)
1226  fall.close()
pydfnworks.dfnGen.meshing.mesh_dfn_helper
Definition: mesh_dfn_helper.py:1
pydfnworks.dfnGen.well_package.wells.find_segments
def find_segments(well)
Definition: wells.py:646
pydfnworks.dfnGen.well_package.wells.well_point_of_intersection
def well_point_of_intersection(well)
Definition: wells.py:800
pydfnworks.dfnGen.well_package.wells.rotation_matrix
def rotation_matrix(normalA, normalB)
Definition: wells.py:1029
pydfnworks.dfnGen.well_package.wells.get_center
def get_center(fracture_id)
Definition: wells.py:1003
pydfnworks.dfnGen.well_package.wells.get_segments
def get_segments(well_line_file)
Definition: wells.py:929
pydfnworks.dfnGen.well_package.wells.cross_check_pts
def cross_check_pts(h)
Definition: wells.py:878
pydfnworks.dfnGen.well_package.wells.expand_well
def expand_well(well)
Definition: wells.py:188
pydfnworks.dfnGen.well_package.wells.cleanup_wells
def cleanup_wells(self, wells)
Definition: wells.py:1127
pydfnworks.dfnGen.well_package.wells.get_normal
def get_normal(fracture_id)
Definition: wells.py:981
pydfnworks.dfnGen.well_package.wells.combine_well_boundary_zones
def combine_well_boundary_zones(self, wells)
Definition: wells.py:1181
pydfnworks.dfnGen.well_package.wells.rotate_point
def rotate_point(p, R)
Definition: wells.py:1102
pydfnworks.dfnGen.well_package.wells.find_well_intersection_points
def find_well_intersection_points(self, wells)
Definition: wells.py:544
pydfnworks.dfnGen.well_package.wells.convert_well_to_polyline_avs
def convert_well_to_polyline_avs(well, h=None)
Definition: wells.py:96
pydfnworks.dfnGen.well_package.wells.get_well_zone
def get_well_zone(well, inp_file)
Definition: wells.py:457
pydfnworks.dfnGen.well_package.wells.run_find_well_intersection_points
def run_find_well_intersection_points(well, h)
Definition: wells.py:600
pydfnworks.dfnGen.well_package.wells.tag_well_in_mesh
def tag_well_in_mesh(self, wells)
Definition: wells.py:9