Add pygeosx_tools python package (#1362)

Adding the pygeosx_tools package

Author: cssherman

pygeosx_tools_package/pygeosx_tools/__init__.py

@@ -0,0 +1,83 @@
+
+import os
+import numpy as np
+
+
+def save_tables(axes, properties, table_root='./tables', axes_names=[]):
+    """
+    Saves a set of tables in GEOSX format
+    Notes: The shape of these arrays should match the length of each axis in the specified order
+           The output directory will be created if it does not exist yet
+           If axes_names are not supplied, then they will be selected based
+           on the dimensionality of the grid: 1D=[t]; 3D=[x, y, z]; 4D=[x, y, z, t]
+
+    Args:
+        axes (list): A list of numpy ndarrays defining the table axes
+        properties (dict): A dict of numpy ndarrays defning the table values
+        table_root (str): The root path for the output directory
+        axes_names (list): A list of names for each potential axis (optional)
+    """
+    # Check to see if the axes, properties have consistent shapes
+    axes_size = tuple([len(x) for x in axes])
+    axes_dimension = len(axes_size)
+    for k, p in properties.items():
+        property_size = np.shape(p)
+        if (property_size != axes_size):
+            print('Property:', k)
+            print('Grid size:', axes_size)
+            print('Property size', property_size)
+            raise Exception('Table dimensions do not match proprerties')
+
+    # Check the axes names
+    if axes_names:
+        if (axes_dimension != len(axes_names)):
+            print('Axes dimensions:', axes_dimension)
+            print('Number of axis names provided:', len(axes_names))
+            raise Exception('The grid dimensions and axes names do not match')
+    else:
+        if (axes_dimension == 1):
+            axes_names = ['t']
+        elif (axes_dimension == 3):
+            axes_names = ['x', 'y', 'z']
+        elif (axes_dimension == 4):
+            axes_names = ['x', 'y', 'z', 't']
+        else:
+            axes_names = ['x%i' % (ii) for ii in range(axes_dimension)]
+
+    # Write the axes
+    os.makedirs(table_root, exist_ok=True)
+    for g, a in zip(axes, axes_names):
+        np.savetxt('%s/%s.csv' % (table_root, a),
+                   g,
+                   fmt='%1.5f',
+                   delimiter=',')
+
+    for k, p in properties.items():
+        np.savetxt('%s/%s.csv' % (table_root, k),
+                   np.reshape(p, (-1), order='F'),
+                   fmt='%1.5e',
+                   delimiter=',')
+
+
+def load_tables(axes_names, property_names, table_root='./tables', extension='csv'):
+    """
+    Load a set of tables in GEOSX format
+
+    Args:
+        axes_names (list): Axis file names in the target directory (with no extension)
+        property_names (list): Property file names in the target directory (with not extension)
+        table_root (str): Root path for the table directory
+        extension (str): Table file extension (default = 'csv')
+
+    Returns:
+        tuple: List of axes values, and dictionary of table values
+    """
+    # Load axes
+    axes = [np.loadtxt('%s/%s.%s' % (table_root, axis, extension), unpack=True) for axis in axes_names]
+    N = tuple([len(x) for x in axes])
+
+    # Load properties
+    properties = {p: np.reshape(np.loadtxt('%s/%s.%s' % (table_root, p, extension)), N, order='F') for p in property_names}
+
+    return axes, properties
+

pygeosx_tools_package/pygeosx_tools/file_io.py

@@ -0,0 +1,108 @@
+
+import numpy as np
+from scipy import stats
+
+
+def apply_to_bins(fn, position, value, bins, collapse_edges=True):
+    """
+    Apply a function to values that are located within a series of bins
+    Note: if a bin is empty, this function will fill a nan value
+
+    Args:
+        fn (function): Function that takes a single scalar or array input
+        position (np.ndarray): A 1D list/array describing the location of each sample
+        value (np.ndarray): A 1D list/array of values at each location
+        bins (np.ndarray): The bin edges for the position data
+        collapse_edges (bool): Controls the behavior of edge-data (default=True)
+
+    Returns:
+        np.ndarray: an array of function results for each bin
+    """
+    # Sort values into bins
+    Nr = len(bins) + 1
+    Ibin = np.digitize(position, bins)
+
+    if collapse_edges:
+        Nr -= 2
+        Ibin -= 1
+        Ibin[Ibin == -1] = 0
+        Ibin[Ibin == Nr] = Nr-1
+
+    # Apply functions to bins
+    binned_values = np.zeros(Nr)
+    for ii in range(Nr):
+        tmp = (Ibin == ii)
+        if np.sum(tmp):
+            binned_values[ii] = fn(value[tmp])
+        else:
+            # Empty bin
+            binned_values[ii] = np.NaN
+
+    return binned_values
+
+
+def extrapolate_nan_values(x, y, slope_scale=0.0):
+    """
+    Fill in any nan values in two 1D arrays by extrapolating
+
+    Args:
+        x (np.ndarray): 1D list/array of positions
+        y (np.ndarray): 1D list/array of values
+        slope_scale (float): value to scale the extrapolation slope (default=0.0)
+
+    Returns:
+        np.ndarray: The input array with nan values replaced by extrapolated data
+    """
+    Inan = np.isnan(y)
+    reg = stats.linregress(x[~Inan], y[~Inan])
+    y[Inan] = reg[0] * x[Inan] * slope_scale + reg[1]
+    return y
+
+
+def get_random_realization(x, bins, value, rand_fill=0, rand_scale=0, slope_scale=0):
+    """
+    Get a random realization for a noisy signal with a set of bins
+
+    Args:
+        x (np.ndarray): 1D list/array of positions
+        bins (np.ndarray): 1D list/array of bin edges
+        value (np.ndarray): 1D list/array of values
+        rand_fill (float): The standard deviation to use where data is not defined (default=0)
+        rand_scale (float): Value to scale the standard deviation for the realization (default=0)
+        slope_scale (float): Value to scale the extrapolation slope (default=0.0)
+
+    Returns:
+        np.ndarray: An array containing the random realization
+    """
+    y_mean = apply_to_bins(np.mean, x, value, bins)
+    y_std = apply_to_bins(np.std, x, value, bins)
+
+    # Extrapolate to fill the upper/lower bounds
+    x_mid = bins[:-1] + 0.5 * (bins[1] - bins[0])
+    y_mean = extrapolate_nan_values(x_mid, y_mean, slope_scale)
+    y_std[np.isnan(y_std)] = rand_fill
+
+    # Add a random perturbation to the target value to match missing high/lows
+    y_final = y_mean + (rand_scale * y_std * np.random.randn(len(y_mean)))
+    return y_final
+
+
+def get_realizations(x, bins, targets):
+    """
+    Get random realizations for noisy signals on target bins
+
+    Args:
+        x (np.ndarray): 1D list/array of positions
+        bins (np.ndarray): 1D list/array of bin edges
+        targets (dict): Dict of geosx target keys, inputs to get_random_realization
+
+    Returns:
+        dict: Dictionary of random realizations
+
+    """
+    results = {}
+    for k, t in targets.items():
+        results[k] = get_random_realization(x, bins, **t)
+    return results
+
+

pygeosx_tools_package/pygeosx_tools/mesh_interpolation.py

@@ -0,0 +1,99 @@
+
+import numpy as np
+import re
+
+
+def parse_las(fname, variable_start='~C', body_start='~A'):
+    """
+    Parse an las format log file
+
+    Args:
+        fname (str): Path to the log file
+        variable_start (str): A string that indicates the start of variable header information (default = '~CURVE INFORMATION')
+        body_start (str): a string that indicates the start of the log body (default = '~A')
+
+    Returns:
+        np.ndarray: a dict containing the values and unit definitions for each variable in the log
+    """
+    results = {}
+    variable_order = []
+
+    # The expected format of the varible definition block is:
+    # name.units code:description
+    variable_regex = re.compile('\s*([^\.^\s]*)\s*(\.[^ ]*) ([^:]*):(.*)')
+
+    with open(fname) as f:
+        file_location = 0
+        for line in f:
+            line = line.split('#')[0]
+            if line:
+                # Preamble
+                if (file_location == 0):
+                    if variable_start in line:
+                        file_location += 1
+
+                # Variable definitions
+                elif (file_location == 1):
+                    # This is not a comment line
+                    if body_start in line:
+                        file_location += 1
+                    else:
+                        match = variable_regex.match(line)
+                        if match:
+                            variable_order.append(match[1])
+                            results[match[1]] = {'units': match[2][0:],
+                                                 'code': match[3],
+                                                 'description': match[4],
+                                                 'values': []}
+                        else:
+                            # As a fall-back use the full line
+                            variable_order.append(line[:-1])
+                            results[line[:-1]] = {'units': '',
+                                                  'code': '',
+                                                  'description': '',
+                                                  'values': []}
+
+                # Body
+                else:
+                    for k, v in zip(variable_order, line.split()):
+                        results[k]['values'].append(float(v))
+
+    # Convert values to numpy arrays
+    for k in results:
+        results[k]['values'] = np.array(results[k]['values'])
+
+    return results
+
+
+def convert_E_nu_to_K_G(E, nu):
+    """
+    Convert young's modulus and poisson's ratio to bulk and shear modulus
+
+    Args:
+        E (float, np.ndarray): Young's modulus
+        nu (float, np.ndarray): Poisson's ratio
+
+    Returns:
+        tuple: bulk modulus, shear modulus with same size as inputs
+    """
+    K = E / (3.0 * (1 - 2.0 * nu))
+    G = E / (2.0 * (1 + nu))
+    return K, G
+
+
+def estimate_shmin(z, rho, nu):
+    """
+    Estimate the minimum horizontal stress using the poisson's ratio
+
+    Args:
+        z (float, np.ndarray): Depth
+        rho (float, np.ndarray): Density
+        nu (float, np.ndarray): Poisson's ratio
+
+    Returns:
+        float: minimum horizontal stress
+    """
+    k = nu / (1.0 - nu)
+    sigma_h = k * rho * 9.81 * z
+    return sigma_h
+