#!/usr/bin/env python
# Calculate mass balance across a boundary of a DFN mesh
# Make sure that you run PFLOTRAN with MASS_FLOWRATE under OUTPUT
# Satish Karra
# March 17, 2016
# Updated Jeffrey Hyman Dec 19 2018
import numpy as np
import glob
from pydfnworks.dfnGen.meshing.mesh_dfn_helper import parse_params_file
__author__ = 'Satish Karra'
__email__ = 'satkarra@lanl.gov'
def get_domain():
''' Return dictionary of domain x,y,z by calling parse_params_file
Parameters
----------
None
Returns
-------
domain : dict
Dictionary of domain sizes in x, y, z
Notes
-----
parse_params_file() is in mesh_dfn_helper.py
'''
_, _, _, _, domain = parse_params_file(quiet=True)
return domain
def parse_pflotran_input(pflotran_input_file):
''' Walk through PFLOTRAN input file and find inflow boundary, inflow and outflow pressure, and direction of flow
Parameters
----------
pflotran_input_file : string
Name of PFLOTRAN input file
Returns
-------
inflow_pressure : double
Inflow Pressure boundary condition
outflow_pressure : float
Outflow pressure boundary condition
inflow_file : string
Name of inflow boundary .ex file
direction : string
Primary direction of flow x, y, or z
Notes
-----
Currently only works for Dirichlet Boundary Conditions
'''
with open(pflotran_input_file) as fp:
outflow_found = False
inflow_found = False
for line in fp.readlines():
if "BOUNDARY_CONDITION OUTFLOW" in line:
outflow_found = True
if outflow_found:
if "REGION" in line:
outflow = line.split()[-1]
outflow_found = False
if "BOUNDARY_CONDITION INFLOW" in line:
inflow_found = True
if inflow_found:
if "REGION" in line:
inflow = line.split()[-1]
inflow_found = False
with open(pflotran_input_file) as fp:
inflow_name_found = False
outflow_name_found = False
inflow_found = False
outflow_found = False
for line in fp.readlines():
if "REGION " + inflow in line:
inflow_name_found = True
if inflow_name_found:
if "FILE" in line:
inflow_file = line.split()[-1]
inflow_name_found = False
if "FLOW_CONDITION " + inflow in line:
inflow_found = True
if inflow_found:
if "PRESSURE " in line:
if "dirichlet" not in line:
tmp = line.split()[-1]
tmp = tmp.split('d')
inflow_pressure = float(tmp[0]) * 10**float(tmp[1])
inflow_found = False
if "REGION " + outflow in line:
outflow_name_found = True
if outflow_name_found:
if "FILE" in line:
outflow_file = line.split()[-1]
outflow_name_found = False
if "FLOW_CONDITION " + outflow in line:
outflow_found = True
if outflow_found:
if "PRESSURE " in line:
if "dirichlet" not in line:
tmp = line.split()[-1]
tmp = tmp.split('d')
outflow_pressure = float(tmp[0]) * 10**float(tmp[1])
outflow_found = False
if inflow_file == 'pboundary_left_w.ex' or inflow_file == 'pboundary_right_e.ex':
direction = 'x'
if inflow_file == 'pboundary_front_n.ex' or inflow_file == 'pboundary_back_s.ex':
direction = 'y'
if inflow_file == 'pboundary_top.ex' or inflow_file == 'pboundary_bottom.ex':
direction = 'z'
print("Inflow file: %s" % inflow_file)
print("Inflow Pressure %e" % inflow_pressure)
print("Outflow Pressure %e" % outflow_pressure)
print("Primary Flow Direction : %s" % direction)
return inflow_pressure, outflow_pressure, inflow_file, direction
def flow_rate(darcy_vel_file, boundary_file):
'''Calculates the flow rate across the inflow boundary
Parameters
----------
darcy_vel_file : string
Name of concatenated Darcy velocity file
boundary_file : string
ex file for the inflow boundary
Returns
-------
mass_rate : float
Mass flow rate across the inflow boundary
volume_rate : float
Volumetric flow rate across the inflow boundary
Notes
-----
None
'''
# Calculate the mass flow rate
mass_rate = 0.0 #kg/s
volume_rate = 0.0 #m^3/s
dat_boundary = np.genfromtxt(boundary_file, skip_header=1)
dat = np.genfromtxt(darcy_vel_file)
for cell in dat_boundary[:, 0]:
if (np.any(dat[:, 0] == int(cell))):
ids = np.where(dat[:, 0] == int(cell))[0]
for idx in ids:
cell_up = int(dat[idx, 0])
cell_down = int(dat[idx, 1])
mass_flux = dat[idx, 2] # in m/s , darcy flux, right? m3/m2/s
density = dat[idx, 3] # in kg/m3
area = dat[idx, 4] # in m^2
if (cell_up == int(cell)):
mass_rate = mass_flux * area * \
density + mass_rate # in kg/s
volume_rate = mass_flux * area + volume_rate #in m3/s
else:
mass_rate = - mass_flux * area * \
density + mass_rate # in kg/s
volume_rate = -mass_flux * area + volume_rate #in m3/s
#print cell_up, cell_down, mass_flux, density, area, mass_rate, volume_rate
if (np.any(dat[:, 1] == int(cell))):
ids = np.where(dat[:, 1] == int(cell))[0]
for idx in ids:
cell_up = int(dat[idx, 0])
cell_down = int(dat[idx, 1])
mass_flux = dat[idx, 2] # in m/s
density = dat[idx, 3] # in kg/m3
area = dat[idx, 4] # in m^2
if (cell_up == int(cell)):
mass_rate = mass_flux * area * \
density + mass_rate # in kg/s
volume_rate = mass_flux * area + volume_rate #in m3/s
else:
mass_rate = - mass_flux * area * \
density + mass_rate # in kg/s
volume_rate = -mass_flux * area + volume_rate #in m3/s
#print cell_up, cell_down, mass_flux, density, area, mass_rate, volume_rate
return mass_rate, volume_rate
def dump_effective_perm(local_jobname, mass_rate, volume_rate, domain,
direction, inflow_pressure, outflow_pressure):
'''Compute the effective permeability of the DFN and write it to screen and to the file local_jobname_effective_perm.dat
Parameters
----------
local_jobname : string
Jobname
mass_rate : float
Mass flow rate through inflow boundary
volume_rate : float
Volumetric flow rate through inflow boundary
direction : string
Primary direction of flow (x, y, or z)
domain : dict
Dictionary of domain sizes in x, y, z
inflow_pressure : float
Inflow boundary pressure
outflow_pressure : float
Outflow boundary pressure
Returns
-------
None
Notes
-----
Information is written into (local_jobname)_effective_perm.txt
'''
Lm3 = domain['x'] * domain['y'] * domain['z'] #L/m^3
# if flow is in x direction, cross section is domain['y']*domain['z']
if direction == 'x':
surface = domain['y'] * domain['z']
L = domain['x']
if direction == 'y':
surface = domain['x'] * domain['z']
L = domain['y']
if direction == 'z':
surface = domain['x'] * domain['y']
L = domain['z']
# Print the calculated mass flow rate in kg/s
mu = 8.9e-4 #dynamic viscosity of water at 20 degrees C, Pa*s
spery = 3600. * 24. * 365.25 #seconds per year
# compute pressure gradient
pgrad = (inflow_pressure - outflow_pressure) / L
print("\n\nEffective Permeabilty Properties\n")
fp = open(local_jobname + '_effective_perm.txt', "w")
print('The mass flow rate [kg/s]: ' + str(mass_rate))
fp.write('The mass flow rate [kg/s]: %e\n' % (mass_rate))
print('The volume flow rate [m3/s]: ' + str(volume_rate))
fp.write('The volume flow rate [m3/s]: %s\n' % (volume_rate))
q = volume_rate / surface #darcy flux over entire face m3/m2/s
if direction == 'x':
print('The darcy flow rate over %f x %f m2 area [m3/m2/s]: %e' %
(domain['y'], domain['z'], q))
fp.write('The darcy flow rate over %f x %f m2 area [m3/m2/s]: %e\n' %
(domain['y'], domain['z'], q))
print('The darcy flow rate over %f x %f m2 area [m3/m2/y]: %e' %
(domain['y'], domain['z'], spery * q))
fp.write('The darcy flow rate over %f x %f m2 area [m3/m2/y]: %e\n' %
(domain['y'], domain['z'], spery * q))
if direction == 'y':
print('The darcy flow rate over %f x %f m2 area [m3/m2/s]: %e' %
(domain['x'], domain['z'], q))
fp.write('The darcy flow rate over %f x %f m2 area [m3/m2/s]: %e\n' %
(domain['x'], domain['z'], q))
print('The darcy flow rate over %f x %f m2 area [m3/m2/y]: %e' %
(domain['x'], domain['z'], spery * q))
fp.write('The darcy flow rate over %f x %f m2 area [m3/m2/y]: %e\n' %
(domain['x'], domain['z'], spery * q))
if direction == 'z':
print('The darcy flow rate over %f x %f m2 area [m3/m2/s]: %e' %
(domain['x'], domain['y'], q))
fp.write('The darcy flow rate over %f x %f m2 area [m3/m2/s]: %e\n' %
(domain['x'], domain['y'], q))
print('The darcy flow rate over %f x %f m2 area [m3/m2/y]: %e' %
(domain['x'], domain['y'], spery * q))
fp.write('The darcy flow rate over %f x %f m2 area [m3/m2/y]: %e\n' %
(domain['x'], domain['y'], spery * q))
print('The effective permeability of the domain [m2]: ' +
str(q * mu / pgrad))
fp.write('The effective permeability of the domain [m2]: %e\n' %
(q * mu / pgrad))
fp.close()
[docs]def effective_perm(self):
'''Computes the effective permeability of a DFN in the primary direction of flow using a steady-state PFLOTRAN solution.
Parameters
----------
self : object
DFN Class
Returns
-------
None
Notes
-----
1. Information is written to screen and to the file self.local_jobname_effective_perm.txt
2. Currently, only PFLOTRAN solutions are supported
3. Assumes density of water
'''
print("\n--> Computing Effective Permeability of Block")
if not self.flow_solver == "PFLOTRAN":
print(
"Incorrect flow solver selected. Cannot compute effective permeability"
)
return 0
darcy_vel_file = 'darcyvel.dat'
pflotran_input_file = self.local_dfnFlow_file
inflow_pressure, outflow_pressure, boundary_file, direction = parse_pflotran_input(
pflotran_input_file)
domain = get_domain()
mass_rate, volume_rate = flow_rate(darcy_vel_file, boundary_file)
dump_effective_perm(self.local_jobname, mass_rate, volume_rate, domain,
direction, inflow_pressure, outflow_pressure)
print("\n--> Complete\n\n")