Source code for pydfnworks.dfnGraph.graph_flow

import networkx as nx
import numpy as np
import sys
import scipy.sparse
import h5py

# pydfnworks modules
from pydfnworks.dfnGraph.intersection_graph import create_intersection_graph
from pydfnworks.dfnGraph.graph_attributes import add_perm, add_area, add_weight


def get_laplacian_sparse_mat(G,
                             nodelist=None,
                             weight=None,
                             dtype=None,
                             format='lil'):
    """ Get the matrices D, A that make up the Laplacian sparse matrix in desired sparsity format. Used to enforce boundary conditions by modifying rows of L = D - A

    Parameters
    ----------
        G : object
            NetworkX graph equipped with weight attribute

        nodelist : list
            list of nodes of G for which laplacian is desired. Default is None in which case, all the nodes
        
        weight : string
            For weighted Laplacian, else all weights assumed unity
        
        dtype :  default is None, cooresponds to float
        
        format: string
            sparse matrix format, csr, csc, coo, lil_matrix with default being lil

    Returns
    -------
        D : sparse 2d float array       
            Diagonal part of Laplacian
            
        A : sparse 2d float array
            Adjacency matrix of graph
    """

    A = nx.to_scipy_sparse_array(G,
                                  nodelist=nodelist,
                                  weight=weight,
                                  dtype=dtype,
                                  format=format)

    (n, n) = A.shape
    data = np.asarray(A.sum(axis=1).T)
    D = scipy.sparse.spdiags(data, 0, n, n, format=format)
    return D, A


def prepare_graph_with_attributes(inflow, outflow, G=None):
    """ Create a NetworkX graph, prepare it for flow solve by equipping edges with  attributes, renumber vertices, and tag vertices which are on inlet or outlet
    
    Parameters
    ----------
        inflow : string
            name of file containing list of DFN fractures on inflow boundary

        outflow: string
            name of file containing list of DFN fractures on outflow boundary

    Returns
    -------
        Gtilde : NetworkX graph
    """

    if G == None:
        G = create_intersection_graph(inflow, outflow)

        Gtilde = G.copy()
        # need to add aperture
        add_perm(Gtilde)
        add_area(Gtilde)
        add_weight(Gtilde)

    else:
        Gtilde = G
        #add_perm(Gtilde)
        #add_area(Gtilde)
        add_weight(Gtilde)

    for v in nx.nodes(Gtilde):
        Gtilde.nodes[v]['inletflag'] = False
        Gtilde.nodes[v]['outletflag'] = False

    if len(list(nx.neighbors(Gtilde, 's'))) == 0:
        error = "Error. There are no nodes in the inlet.\nExiting"
        sys.stderr.write(error)
        sys.exit(1)

    for v in nx.neighbors(Gtilde, 's'):
        Gtilde.nodes[v]['inletflag'] = True

    if len(list(nx.neighbors(Gtilde, 't'))) == 0:
        error = "Error. There are no nodes in the outlet.\nExiting"
        sys.stderr.write(error)
        sys.exit(1)

    for v in nx.neighbors(Gtilde, 't'):
        Gtilde.nodes[v]['outletflag'] = True

    Gtilde.remove_node('s')
    Gtilde.remove_node('t')

    Gtilde = nx.convert_node_labels_to_integers(Gtilde,
                                                first_label=0,
                                                ordering="sorted",
                                                label_attribute="old_label")

    return Gtilde


def solve_flow_on_graph(G, pressure_in, pressure_out, fluid_viscosity, phi):
    """ Given a NetworkX graph prepared  for flow solve, solve for vertex pressures, and equip edges with attributes (Darcy) flux  and time of travel

    Parameters
    ----------
        G : NetworkX graph

        pressure_in : double
            Value of pressure (in Pa) at inlet
        
        pressure_out : double
            Value of pressure (in Pa) at outlet
        
        fluid_viscosity : double
            optional, in Pa-s, default is for water

        phi : double
            Porosity, default is 1

    Returns
    -------
        H : Acyclic Directed NetworkX graph 
            H is updated with vertex pressures, edge fluxes and travel times. The only edges that exists are those with postive flow rates. 

    Notes
    ----------
        None

    """

    print("--> Starting Graph flow")

    Inlet = [v for v in nx.nodes(G) if G.nodes[v]['inletflag']]
    Outlet = [v for v in nx.nodes(G) if G.nodes[v]['outletflag']]

    if not set(Inlet).isdisjoint(set(Outlet)):
        error = "Incompatible graph: Vertex connected to both source and target\n"
        sys.stderr.write(error)
        sys.exit(1)

    D, A = get_laplacian_sparse_mat(G, weight='weight', format='lil')

    b = np.zeros(G.number_of_nodes())

    for v in Inlet:
        b[v] = pressure_in
        A[v, :] = 0
        D[v, v] = 1.0
    for v in Outlet:
        b[v] = pressure_out
        A[v, :] = 0
        D[v, v] = 1.0
    L = D - A  # automatically converts to csr when returning L

    print("--> Solving Linear System for pressure at nodes")
    pressure = scipy.sparse.linalg.spsolve(L, b)
    print("--> Updating graph edges with flow solution")

    for v in nx.nodes(G):
        G.nodes[v]['pressure'] = pressure[v]

    H = nx.DiGraph()
    H.add_nodes_from(G.nodes(data=True))

    for u, v in nx.edges(G):
        # Find direction of flow
        if G.nodes[u]['pressure'] > G.nodes[v]['pressure']:
            upstream = u
            downstream = v
        elif G.nodes[v]['pressure'] >= G.nodes[u]['pressure']:
            upstream = v
            downstream = u

        delta_p = G.nodes[upstream]['pressure'] - G.nodes[downstream][
            'pressure']
        if delta_p > 1e-16:
            ## Create new edge in DiGraph
            H.add_edge(upstream, downstream)
            # Transfer edge attributes
            for att in [
                    'perm', 'iperm', 'length', 'weight', 'area', 'frac', 'b'
            ]:
                H.edges[upstream, downstream][att] = G.edges[upstream,
                                                             downstream][att]

            H.edges[upstream, downstream]['flux'] = (
                H.edges[upstream, downstream]['perm'] / fluid_viscosity) * (
                    delta_p / H.edges[upstream, downstream]['length'])

            H.edges[upstream, downstream]['vol_flow_rate'] = H.edges[
                upstream, downstream]['flux'] * H.edges[upstream,
                                                        downstream]['area']

            # H.edges[downstream, upstream]['vol_flow_rate'] =  -1*H.edges[upstream, downstream]['vol_flow_rate']

            H.edges[upstream,
                    downstream]['velocity'] = H.edges[upstream,
                                                      downstream]['flux'] / phi
            H.edges[upstream,
                    downstream]['time'] = H.edges[upstream, downstream][
                        'length'] / (H.edges[upstream, downstream]['velocity'])

    print("--> Graph flow complete")
    return H


def compute_dQ(self, G):
    """ Computes the DFN fracture intensity (p32) and flow channeling density indicator from the graph flow solution on G

    Parameters
    -----------------
        self : object
            DFN Class

        G : networkX graph 
            Output of run_graph_flow

    Returns
    ---------------
        p32 : float
            Fracture intensity
        dQ : float flow channeling density indicator 

    Notes
    ------------
        For definitions of p32 and dQ along with a discussion see " Hyman, Jeffrey D. "Flow channeling in fracture networks: characterizing the effect of density on preferential flow path formation." Water Resources Research 56.9 (2020): e2020WR027986. "

    """
    print(
        "--> Computing fracture intensity (p32) and flow channeling density indicator (dQ)"
    )

    fracture_surface_area = 2*self.surface_area
    domain_volume = self.domain['x'] * self.domain['y'] * self.domain['z']

    Qf = np.zeros(self.num_frac)
    ## convert to undirected
    H = G.to_undirected()
    ## walk through fractures
    for curr_frac in range(1, self.num_frac + 1):
        # print(f"\nstarting on fracture {curr_frac}")
        # Gather nodes on current fracture
        current_nodes = []
        for u, d in H.nodes(data=True):
            for f in d["frac"]:
                if f == curr_frac:
                    current_nodes.append(u)
        # cycle through nodes on the fracture and get the outgoing / incoming
        # volumetric flow rates
        for u in current_nodes:
            neighbors = H.neighbors(u)
            for v in neighbors:
                if v not in current_nodes:
                    # outgoing vol flow rate
                    Qf[curr_frac - 1] += abs(H[u][v]['vol_flow_rate'])
                    for f in H.nodes[v]['frac']:
                        if f != curr_frac and f != 's' and f != 't':
                            # incoming vol flow rate
                            Qf[f - 1] += abs(H[u][v]['vol_flow_rate'])
    # Divide by 1/2 to remove up double counting
    Qf *= 0.5
    p32 = fracture_surface_area.sum() / domain_volume
    top = sum(fracture_surface_area * Qf)**2
    bottom = sum(fracture_surface_area * Qf**2)
    print(top, bottom)
    dQ = (1.0 / domain_volume) * (top / bottom)
    print(f"--> P32: {p32:0.2e} [1/m]")
    print(f"--> dQ: {dQ:0.2e} [1/m]")
    print(f"--> Active surface percentage {100*dQ/p32:0.2f}")
    print(f"--> Geometric equivalent fracture spacing {1/p32:0.2e} m")
    print(f"--> Hydrological equivalent fracture spacing {1/dQ:0.2e} m")
    print("--> Complete \n")
    return p32, dQ, Qf


def dump_graph_flow_values(G):

    num_edges = G.number_of_edges()
    velocity = np.zeros(num_edges)
    lengths = np.zeros_like(velocity)
    vol_flow_rate = np.zeros_like(velocity)
    area = np.zeros_like(velocity)
    aperture = np.zeros_like(velocity)
    volume = np.zeros_like(velocity)

    for i, val in enumerate(G.edges(data=True)):
        u, v, d = val
        velocity[i] = d['velocity']
        lengths[i] = d['length']
        vol_flow_rate[i] = d['vol_flow_rate']
        area[i] = d['area']
        aperture[i] = d['b']
        volume[i] = area[i] * aperture[i]

    with h5py.File(f"graph_flow.hdf5", "w") as f5file:
        h5dset = f5file.create_dataset('velocity', data=velocity)
        h5dset = f5file.create_dataset('length', data=lengths)
        h5dset = f5file.create_dataset('vol_flow_rate', data=vol_flow_rate)
        h5dset = f5file.create_dataset('area', data=area)
        h5dset = f5file.create_dataset('aperture', data=aperture)
        h5dset = f5file.create_dataset('volume', data=volume)
    f5file.close()


[docs] def run_graph_flow(self, inflow, outflow, pressure_in, pressure_out, fluid_viscosity=8.9e-4, phi=1, G=None): """ Solve for pressure driven steady state flow on a graph representation of the DFN. Parameters ---------- self : object DFN Class inflow : string name of file containing list of DFN fractures on inflow boundary outflow: string name of file containing list of DFN fractures on outflow boundary pressure_in : double Value of pressure at inlet [Pa] pressure_out : double Value of pressure at outlet [Pa] fluid_viscosity : double optional, default is for water. [Pa*s] phi : double Fracture porosity, default is 1 [-] G : Input Graph Returns ------- Gtilde : NetworkX graph Gtilde is a directed acyclic graph with vertex pressures, fluxes, velocities, volumetric flow rates, and travel times """ if G == None: G = self.create_graph("intersection", inflow, outflow) Gtilde = prepare_graph_with_attributes(inflow, outflow, G) Gtilde = solve_flow_on_graph(Gtilde, pressure_in, pressure_out, fluid_viscosity, phi) dump_graph_flow_values(Gtilde) return Gtilde