upscale.py

Go to the documentation of this file.

"""
.. module:: upscale.py
   :synopsis: Generates PFLOTRAN or FEHM input from octree refined continuum mesh
.. moduleauthor:: Matthew Sweeney <msweeney2796@lanl.gov>
 
"""
 
import os
import numpy as np
import sys
import subprocess
import shutil
import h5py
from pydfnworks.dfnGen.meshing import mesh_dfn_helper as mh
import time
import math as m
import glob
import pickle
 
 
def upscale(self, mat_perm, mat_por, path='../'):
    """ Generate permeabilities and porosities based on output of map2continuum.
 
    Parameters
    ----------
        self : object
            DFN Class 
        mat_perm : float 
            Matrix permeability (in m^2)
        mat_por: float
            Matrix porosity
 
    Returns
    -------
        perm_fehm.dat : text file
            Contains permeability data for FEHM input
        rock_fehm.dat : text file
            Contains rock properties data for FEHM input
        mesh_permeability.h5 : h5 file
            Contains permeabilites at each node for PFLOTRAN input
        mesh_porosity.h5 : h5 file
            Contains porosities at each node for PFLOTRAN input               
        
    Notes
    -----
        None
 
    """
 
    print('=' * 80)
    print("Generating permeability and porosity for octree mesh: Starting")
    print('=' * 80)
 
    # Bring in all the relevant data
    try:
        os.symlink(path + 'params.txt', 'params.txt')
    except:
        pass
    num_poly, _, _, _, domain = mh.parse_params_file(quiet=True)
 
    aperture = np.genfromtxt(path + 'aperture.dat', skip_header=1)[:, -1]
    normal_vectors = np.genfromtxt(path + 'normal_vectors.dat', delimiter=' ')
 
    if self.flow_solver == "FEHM":
        with open("perm_fehm.dat", "w") as f:
            f.write("perm\n")
        with open("rock_fehm.dat", "w") as g:
            g.write("rock\n")
 
    # Make dictionary w/ cell IDs (keys) and intersecting fractures, areas  (values)
    for i in range(1, num_poly + 1):
        if i == 1:
            with open('frac1.inp', 'r') as f:
                line = f.readline().strip().split()
                num_nodes = int(float(line[0]))
                num_cells = int(float(line[1]))
                f_dict = {}
                perm_var = np.zeros(num_nodes, 'float')
                por_var = np.zeros(num_nodes, 'float')
                cv_vol = np.zeros(num_nodes, 'float')
                iarray = np.zeros(num_nodes, '=i4')
                frac_vol = np.zeros(num_nodes, 'float')
                permX = np.zeros(num_nodes, 'float')
                permY = np.zeros(num_nodes, 'float')
                permZ = np.zeros(num_nodes, 'float')
            f.close()
        with open('frac{0}.inp'.format(i), 'r') as f:
            with open('area_sum{0}.inp'.format(i), 'r') as g:
                for j in range(num_nodes):
                    f.readline()
                    g.readline()
                f.readline()
                g.readline()
                for j in range(num_cells):
                    f.readline()
                    g.readline()
                f.readline()
                f.readline()
                f.readline()
                g.readline()
                g.readline()
                g.readline()
                g.readline()
                g.readline()
                g.readline()
                g.readline()
                for j in range(num_nodes):
                    fline = f.readline().strip().split()
                    gline = g.readline().strip().split()
                    iarray[j] = int(float(gline[0]))
                    if int(float(gline[1])) != (num_poly + 1) and float(
                            gline[6]) > 0:
                        f_dict.setdefault(j + 1, []).append(
                            (i, float(gline[6])))
            g.close()
        f.close()
 
    p_out = open("connections.p", "wb")
    pickle.dump(f_dict, p_out)
    p_out.close()
 
    with open('full_mesh.uge', 'r') as f:
        f.readline()
        for j in range(num_nodes):
            fline = f.readline().strip().split()
            cv_vol[j] = float(fline[4])
    f.close()
 
    # Populate permeability and porosity arrays here
    for i in range(1, num_nodes + 1):
        if i in f_dict:
 
            # Get porosity:
            for j in range(len(f_dict[i])):
                # Calculate total volume of fractures in cv cell i
                frac_vol[i -
                         1] += aperture[f_dict[i][j][0] - 1] * f_dict[i][j][1]
            por_var[i - 1] = frac_vol[i - 1] / cv_vol[i - 1]
            if por_var[i - 1] == 0:
                por_var[i - 1] = mat_por
            if por_var[i - 1] > 1.0:
                por_var[i - 1] = 1.0
 
            # Get permeability:
            perm_tensor = np.zeros([3, 3])
            phi_sum = 0
            for j in range(len(f_dict[i])):
                phi = (aperture[f_dict[i][j][0] - 1] *
                       f_dict[i][j][1]) / cv_vol[i - 1]
                if phi > 1.0:
                    phi = 1.0
                phi_sum += phi
                if phi_sum > 1.0:
                    phi_sum = 1.0
                b = aperture[f_dict[i][j][0] - 1]
                # Construct tensor Omega
                Omega = np.zeros([3, 3])
                n1 = normal_vectors[f_dict[i][j][0] - 1][0]
                n2 = normal_vectors[f_dict[i][j][0] - 1][1]
                n3 = normal_vectors[f_dict[i][j][0] - 1][2]
                Omega[0][0] = (n2)**2 + (n3)**2
                Omega[0][1] = -n1 * n2
                Omega[0][2] = -n3 * n1
                Omega[1][0] = -n1 * n2
                Omega[1][1] = (n3)**2 + (n1)**2
                Omega[1][2] = -n2 * n3
                Omega[2][0] = -n3 * n1
                Omega[2][1] = -n2 * n3
                Omega[2][2] = (n1)**2 + (n2)**2
                perm_tensor += (phi * (b)**2 * Omega)
            perm_tensor *= 1. / 12
 
            # Calculate eigenvalues
            permX[i - 1] = np.linalg.eigvals(perm_tensor)[0]
            permY[i - 1] = np.linalg.eigvals(perm_tensor)[1]
            permZ[i - 1] = np.linalg.eigvals(perm_tensor)[2]
 
            # Arithmetic average of matrix perm
            permX[i - 1] += (1 - phi_sum) * mat_perm
            permY[i - 1] += (1 - phi_sum) * mat_perm
            permZ[i - 1] += (1 - phi_sum) * mat_perm
 
            # Correction factor
 
            # Actual value doesn't matter here, just needs to be high
            min_n1 = 1e6
            min_n2 = 1e6
            min_n3 = 1e6
 
            # See Sweeney et al. 2019 Computational Geoscience
            for j in range(len(f_dict[i])):
                n1_temp = normal_vectors[f_dict[i][j][0] - 1][0]
                theta1_t = m.degrees(m.acos(n1_temp)) % 90
                if abs(theta1_t - 45) <= min_n1:
                    theta1 = theta1_t
                    min_n1 = theta1_t
                n2_temp = normal_vectors[f_dict[i][j][0] - 1][1]
                theta2_t = m.degrees(m.acos(n2_temp)) % 90
                if abs(theta2_t - 45) <= min_n2:
                    theta2 = theta2_t
                    min_n2 = theta2_t
                n3_temp = normal_vectors[f_dict[i][j][0] - 1][2]
                theta3_t = m.degrees(m.acos(n3_temp)) % 90
                if abs(theta3_t - 45) <= min_n3:
                    theta3 = theta3_t
                    min_n3 = theta3_t
 
            sl = (2 * 2**(1. / 2) - 1) / -45.0
            b = 2 * 2**(1. / 2)
 
            cf_x = sl * abs(theta1 - 45) + b
            cf_y = sl * abs(theta2 - 45) + b
            cf_z = sl * abs(theta3 - 45) + b
 
            permX[i - 1] *= cf_x
            permY[i - 1] *= cf_y
            permZ[i - 1] *= cf_z
 
            # Correct 0 perm if exists
            # Note I don't think it's possible for 0 perm to exist
            # because of L177-9 and 110.  Will leave for now...
            if permX[i - 1] == 0:
                permX[i - 1] += mat_perm
            if permY[i - 1] == 0:
                permY[i - 1] += mat_perm
            if permZ[i - 1] == 0:
                permZ[i - 1] += mat_perm
            perm_var[i - 1] = max(permX[i - 1], permY[i - 1], permZ[i - 1])
        else:
            # Assign matrix properties
            por_var[i - 1] = mat_por
            perm_var[i - 1] = mat_perm
 
            # Note these aren't needed if not using anisotropic perm
            permX[i - 1] = mat_perm
            permY[i - 1] = mat_perm
            permZ[i - 1] = mat_perm
 
        if self.flow_solver == "FEHM":
            with open("perm_fehm.dat", "a") as f:
                f.write(
                    str(i) + " " + str(i) + " " + "1" + " " +
                    str(perm_var[i - 1]) + " " + str(perm_var[i - 1]) + " " +
                    str(perm_var[i - 1]) + "\n")
            with open("rock_fehm.dat", "a") as g:
                g.write(
                    str(i) + " " + str(i) + " " + "1" + " " + "2165." + " " +
                    "931." + " " + str(por_var[i - 1]) + "\n")
 
    # Need an extra space at end for FEHM
    if self.flow_solver == "FEHM":
        with open("perm_fehm.dat", "a") as f:
            f.write("\n")
        with open("rock_fehm.dat", "a") as g:
            g.write("\n")
 
    if self.flow_solver == "PFLOTRAN":
        perm_filename = 'mesh_permeability.h5'
        poros_filename = 'mesh_porosity.h5'
 
        h5file = h5py.File(perm_filename, mode='w')
        dataset_name = 'Cell Ids'
        h5dset = h5file.create_dataset(dataset_name, data=iarray)
 
        dataset_name = 'Permeability'
        h5dset = h5file.create_dataset(dataset_name, data=perm_var)
 
        #dataset_name = 'Perm_X'
        #h5dset = h5file.create_dataset(dataset_name, data = permX)
 
        #dataset_name = 'Perm_Y'
        #h5set = h5file.create_dataset(dataset_name, data = permY)
 
        #dataset_name = 'Perm_Z'
        #h5set = h5file.create_dataset(dataset_name, data = permZ)
        h5file.close()
 
        h5file = h5py.File(poros_filename, mode='w')
        dataset_name = 'Cell Ids'
        h5dset = h5file.create_dataset(dataset_name, data=iarray)
 
        dataset_name = 'Porosity'
        h5dset = h5file.create_dataset(dataset_name, data=por_var)
        h5file.close()
 
    upscale_cleanup()
 
    # What nodes are fractures vs. matrix
    tag = perm_var > mat_perm
    tag.astype("uint8")
    np.savetxt("tag_frac.dat", tag, '%d', delimiter=",")
 
 
def upscale_cleanup():
    files_to_remove = [
        "area*", "build*", "driver*", "ex*", "frac*", "hex*", "intersect*",
        "log*", "out*", "parame*", "remove*", "time*", "tmp*"
    ]
    for name in files_to_remove:
        for fl in glob.glob(name):
            os.remove(fl)