various changes for compatibility with Numpy 1.20

README.md

@@ -38,7 +38,7 @@ Bugs should be reported as issues via GitHub.
 
 ## Getting started
 
-Running PyFrac will require a functioning installation of Python 3, numpy, SciPy and matplotlib. You can download ANACONDA distribution with all these packages from <https://www.anaconda.com/distribution/>. To run PyFrac using Unix shell or windows command prompt, the directory containing the PyFrac source code is required to be added to the ``PYTHONPATH`` environment variable. It can be added with the following command on windows (give the local path of the ``src`` folder in place of path_of_PyFrac_src):
+Running PyFrac will require a functioning installation of Python 3, numpy, SciPy and matplotlib. To run PyFrac using Unix shell or windows command prompt, the directory containing the PyFrac source code is required to be added to the ``PYTHONPATH`` environment variable. It can be added with the following command on windows (give the local path of the ``src`` folder in place of path_of_PyFrac_src):
 
     set PYTHONPATH=%PYTHONPATH%;path_of_PyFrac_src
 

src/HF_reference_solutions.py

@@ -99,8 +99,8 @@ def MDR_M_vertex_solution(Eprime, Q0, density, visc, Mesh, R=None, t=None, requi
             rho[actvElts] ** 3)) / (var1 ** 0.25925925925925924) + (0.07471989686220308 * (rho[actvElts] ** 4)) / (
             var1 ** 0.25925925925925924))
         if np.isinf(p[actvElts]).any():
-            p[p == np.NINF] = min(p[p != np.NINF])
-            p[p == np.Inf] = max(p[p != np.inf])
+            p[p == -np.inf] = min(p[p != -np.inf])
+            p[p == np.inf] = max(p[p != np.inf])
     else:
         p = None
 
@@ -184,8 +184,8 @@ def M_vertex_solution(Eprime, Q0, muPrime, Mesh, inj_point, R=None, t=None, requ
                         0.103558 * rho[actvElts] ** 3 + (1 - rho[actvElts]) ** (1 / 3) * np.log(1 / rho[actvElts]))
                        ) / (t ** (1 / 3) * (1 - rho[actvElts]) ** (1 / 3))
         if np.isinf(p[actvElts]).any():
-            p[p == np.NINF] = min(p[p != np.NINF])
-            p[p == np.Inf] = max(p[p != np.inf])
+            p[p == -np.inf] = min(p[p != -np.inf])
+            p[p == np.inf] = max(p[p != np.inf])
 
     else:
         p = None
@@ -673,8 +673,8 @@ def Mt_vertex_solution(Eprime, Cprime, Q0, muPrime, Mesh, inj_point, R=None, t=N
         Q0 ** 0.125 * var1 * t ** 0.1875)
 
         if np.isinf(p[actvElts]).any():
-            p[p == np.NINF] = min(p[p != np.NINF])
-            p[p == np.Inf] = max(p[p != np.inf])
+            p[p == -np.inf] = min(p[p != -np.inf])
+            p[p == np.inf] = max(p[p != np.inf])
     else:
         p = None
 

src/anisotropy.py

@@ -87,7 +87,7 @@ def find_angle(elt_ribbon, elt_tip, zr_vrtx_tip, a_tip, b_tip, c_tip, x_lft, y_l
     front line in the given tip elements.
     """
 
-    closest_tip_cell = np.zeros((len(elt_ribbon),), dtype=np.int)
+    closest_tip_cell = np.zeros((len(elt_ribbon),), dtype=np.int64)
     dist_ribbon = np.zeros((len(elt_ribbon),), dtype=np.float64)
     alpha = np.zeros((len(elt_ribbon),), dtype=np.float64)
     for i in range(len(elt_ribbon)):
@@ -363,7 +363,7 @@ def construct_polygon(elt_tip, l_tip, alpha_tip, mesh, zero_vertex_tip):
     polygon = np.vstack({tuple(row) for row in polygon})
 
 
-    tip_smoothed = np.array([], dtype=np.int) # the cells containing the edges of polygon (giving the smoothed front)
+    tip_smoothed = np.array([], dtype=np.int64) # the cells containing the edges of polygon (giving the smoothed front)
     smthed_tip_points_left = np.empty((0, 2), dtype=np.float64) #left points of the tip line in the new tip cells
     smthed_tip_points_rgt = np.empty((0, 2), dtype=np.float64) #right points of the tip line in the new tip cells
 
@@ -383,7 +383,7 @@ def construct_polygon(elt_tip, l_tip, alpha_tip, mesh, zero_vertex_tip):
             dist = (polygon[cell_pnt[0], 0] - polygon[cell_pnt, 0]) ** 2 + (polygon[cell_pnt[0], 1] - polygon[
                 cell_pnt, 1]) ** 2
             farthest = np.argmax(dist)
-            to_delete = np.array([], dtype=np.int)
+            to_delete = np.array([], dtype=np.int64)
             for m in range(1, cell_pnt.size):
                 if m != farthest:
                     to_delete = np.append(to_delete, cell_pnt[m])
@@ -421,8 +421,8 @@ def construct_polygon(elt_tip, l_tip, alpha_tip, mesh, zero_vertex_tip):
     smthed_tip_lines_c[zero_angle] = -smthed_tip_points_rgt[zero_angle, 0]
 
     # find the left neighbor of the tip cells in the tip
-    tip_lft_neghb = np.zeros((len(tip_smoothed),), dtype=np.int)
-    tip_rgt_neghb = np.empty((len(tip_smoothed),), dtype=np.int)
+    tip_lft_neghb = np.zeros((len(tip_smoothed),), dtype=np.int64)
+    tip_rgt_neghb = np.empty((len(tip_smoothed),), dtype=np.int64)
     for i in range(len(tip_smoothed)):
         equal = smthed_tip_points_rgt == smthed_tip_points_left[i]
         left_nei = np.where(np.logical_and(equal[:, 0], equal[:, 1]))[0]

src/controller.py

@@ -67,8 +67,8 @@ def __init__(self, Fracture, Solid_prop, Fluid_prop, Injection_prop, Sim_prop, L
         self.stagnant_TS = None     # time step if the front is stagnant. It is increased exponentialy to avoid uneccessary small steps.
         self.perfData = []
         self.lastSavedFile = 0
-        self.lastSavedTime = np.NINF
-        self.lastPlotTime = np.NINF
+        self.lastSavedTime = -np.inf
+        self.lastPlotTime = -np.inf
         self.TmStpCount = 0
         self.chkPntReattmpts = 0    # the number of re-attempts done from the checkpoint. Simulation is declared failed after 5 attempts.
         self.TmStpReductions = 0    # the number of times the time step has been reattempted because the fracture it was advancing too more than two cells in a row

src/mesh.py

@@ -215,7 +215,7 @@ def __init__(self, Lx, Ly, nx, ny, symmetric=False):
         booleconnEdgesNodes = np.zeros([numberofedges, 1], dtype=int)  # Peruzzo 2019
         connEdgesNodes = np.empty([numberofedges, 2], dtype=int)  # Peruzzo 2019
         connElemEdges = np.empty([self.NumberOfElts, 4], dtype=int)  # Peruzzo 2019
-        connEdgesElem = np.full([numberofedges, 2], np.nan, dtype=np.int)  # Peruzzo 2019
+        connEdgesElem = np.full([numberofedges, 2], np.nan, dtype=np.int64)  # Peruzzo 2019
         connNodesEdges = np.full([self.NumberofNodes, 4], np.nan, dtype=int)  # Peruzzo 2019
         connNodesElem = np.full([self.NumberofNodes, 4], np.nan, dtype=int)  # Peruzzo 2019
         k = 0

src/postprocess_fracture.py

@@ -1522,7 +1522,7 @@ def get_velocity_slice(Solid, Fluid, Fr_list, initial_point, vel_direction = 'ux
 
     for i in range(nOFtimes): #each fr has its own mesh
         # 1) get the coordinates of the points in the slices
-        vector_to_be_lost = np.zeros(Fr_list[i].mesh.NumberOfElts,dtype=np.int)
+        vector_to_be_lost = np.zeros(Fr_list[i].mesh.NumberOfElts,dtype=np.int64)
         NotUsd_var_values, sampling_line_center, sampling_cells = get_fracture_variable_slice_cell_center(vector_to_be_lost,
                                                                                                             Fr_list[i].mesh,
                                                                                                             point = initial_point,

src/time_step_solution.py

@@ -987,7 +987,7 @@ def injection_extended_footprint(w_k, Fr_lstTmStp, C, timeStep, Qin, mat_propert
     Fr_kplus1.TarrvlZrVrtx[Fr_kplus1.EltTip[new_tip]] = Fr_kplus1.time - Fr_kplus1.l[new_tip] / Fr_kplus1.v[new_tip]
     Fr_kplus1.wHist = np.maximum(Fr_kplus1.w, Fr_lstTmStp.wHist)
     Fr_kplus1.closed = data[1]
-    tip_neg_rib = np.asarray([], dtype=np.int)
+    tip_neg_rib = np.asarray([], dtype=np.int64)
     # adding tip cells with closed corresponding ribbon cells to the list of closed cells
     for i, elem in enumerate(Fr_kplus1.EltTip):
         if corr_ribbon[i] in Fr_kplus1.closed and elem not in Fr_kplus1.closed:
@@ -1352,7 +1352,7 @@ def solve_width_pressure(Fr_lstTmStp, sim_properties, fluid_properties, mat_prop
             # The code below finds the indices(in the EltCrack list) of the neighbours of all the cells in the crack.
             # This is done to avoid costly slicing of the large numpy arrays while making the linear system during the
             # fixed point iterations. For neighbors that are outside the fracture, len(EltCrack) + 1 is returned.
-            corr_nei = np.full((len(EltCrack_k), 4), len(EltCrack_k), dtype=np.int)
+            corr_nei = np.full((len(EltCrack_k), 4), len(EltCrack_k), dtype=np.int64)
             for i, elem in enumerate(EltCrack_k):
                 corresponding = np.where(EltCrack_k == Fr_lstTmStp.mesh.NeiElements[elem, 0])[0]
                 if len(corresponding) > 0: