diff --git a/t2listing.py b/t2listing.py
index 4c6b8d89..92507246 100755
--- a/t2listing.py
+++ b/t2listing.py
@@ -10,7 +10,6 @@
You should have received a copy of the GNU Lesser General Public License along with PyTOUGH. If not, see ."""
-import string
try:
import numpy as np
from numpy import float64
@@ -52,7 +51,8 @@ def __init__(self, cols, rows, row_format = None, row_line = None, num_keys = 1,
self._row = dict([(r,i) for i,r in enumerate(rows)])
self._data = np.zeros((len(rows), len(cols)), float64)
- def __repr__(self): return repr(self.column_name)+'\n'+repr(self._data)
+ def __repr__(self):
+ return repr(self.column_name) + '\n' + repr(self._data)
def __getitem__(self,key):
if isinstance(key,int):
@@ -92,7 +92,7 @@ def key_from_line(self,line):
if len(key) == 1: return key[0]
else: return tuple(key)
- def rows_matching(self,pattern,index=0,match_any=False):
+ def rows_matching(self, pattern, index = 0, match_any = False):
"""Returns rows in the table with keys matching the specified regular expression pattern
string.
For tables with multiple keys, pattern can be a list or tuple of regular expressions. If
@@ -102,18 +102,20 @@ def rows_matching(self,pattern,index=0,match_any=False):
patterns (instead of all of them). If this option is used in conjunction with a single
string pattern, the specified pattern is applied to all keys."""
from re import search
- if self.num_keys==1: return [self[key] for key in self.row_name if search(pattern,key)]
+ if self.num_keys == 1:
+ return [self[key] for key in self.row_name if search(pattern,key)]
else:
- if isinstance(pattern,str): pattern=[pattern]
- else: pattern=list(pattern)
- if len(pattern)0: keywords.append(self.short_types[0])
- endfile=False
+ self._fullpos, self._pos, self._short = [], [], []
+ fullt, fulls, t, s = [], [], [], []
+ keywords = ['EEEEE']
+ if len(self.short_types) > 0: keywords.append(self.short_types[0])
+ endfile = False
while not endfile:
- kwfound=self.skipto(keywords)
+ kwfound = self.skipto(keywords)
if kwfound:
self._pos.append(self._file.tell())
self.read_header_AUTOUGH2()
- if kwfound=='EEEEE': # full results
+ if kwfound == 'EEEEE': # full results
self._fullpos.append(self._pos[-1])
fullt.append(self.time)
fulls.append(self.step)
@@ -417,11 +427,11 @@ def setup_pos_AUTOUGH2(self):
s.append(self.step)
self._file.readline()
self.skipto(kwfound) # to end of table
- else: endfile=True
- self.times=np.array(t)
- self.steps=np.array(s)
- self.fulltimes=np.array(fullt)
- self.fullsteps=np.array(fulls)
+ else: endfile = True
+ self.times = np.array(t)
+ self.steps = np.array(s)
+ self.fulltimes = np.array(fullt)
+ self.fullsteps = np.array(fulls)
def setup_pos_TOUGH2(self):
"""Sets up _pos list for TOUGH2 listings, containing file position at
@@ -453,44 +463,45 @@ def setup_pos_TOUGH2(self):
def set_table_attributes(self):
"""Makes tables in self._table accessible as attributes."""
- for key,table in self._table.items(): setattr(self,key,table)
+ for key,table in self._table.items(): setattr(self, key, table)
def setup_tables_AUTOUGH2(self):
"""Sets up configuration of element, connection and generation tables."""
- tablename='element'
+ tablename = 'element'
self._file.seek(self._fullpos[0])
while tablename:
self.read_header()
if tablename in self.skip_tables: self.skip_table(tablename)
else: self.setup_table(tablename)
- tablename=self.next_table()
+ tablename = self.next_table()
def setup_tables_TOUGH2(self):
self.read_title()
- tablename='element'
+ tablename = 'element'
self._file.seek(self._fullpos[0])
self.read_header() # only one header at each time
while tablename:
if tablename in self.skip_tables: self.skip_table(tablename)
else: self.setup_table(tablename)
- tablename=self.next_table()
+ tablename = self.next_table()
def setup_tables_TOUGHplus(self):
self.read_title()
- tablename='element'
+ tablename = 'element'
self._file.seek(self._fullpos[0])
self.read_header() # only one header at each time
- nelt_tables=0 # can have multiple element tables
+ nelt_tables = 0 # can have multiple element tables
while tablename:
if tablename in self.skip_tables: self.skip_table(tablename)
else: self.setup_table(tablename)
- tablename=self.next_table()
- if tablename=='element':
- nelt_tables+=1
- tablename+=str(nelt_tables)
+ tablename = self.next_table()
+ if tablename == 'element':
+ nelt_tables += 1
+ tablename += str(nelt_tables)
def next_table_AUTOUGH2(self):
- """Goes to start of next table at current time and returns its type- or None if there are no more."""
+ """Goes to start of next table at current time and returns its type,
+ or None if there are no more."""
keyword = self.readline()[1:6]
return self.table_type(keyword)
@@ -532,7 +543,8 @@ def next_table_TOUGHplus(self):
def skip_to_table_AUTOUGH2(self,tablename,last_tablename,nelt_tables):
"""Skips forwards to headers of table with specified name at the current time."""
tablechar = tablename[0].upper()
- if self._short[self._index]: keyword, first_tablechar = tablechar+'SHORT', self.short_types[0][0]
+ if self._short[self._index]:
+ keyword, first_tablechar = tablechar + 'SHORT', self.short_types[0][0]
else: keyword, first_tablechar = tablechar*5, 'E'
if tablechar != first_tablechar: self.skipto(keyword)
self.skipto('OUTPUT')
@@ -544,31 +556,31 @@ def skip_to_table_TOUGH2(self,tablename,last_tablename,nelt_tables):
if last_tablename is None:
self.skipto('@@@@@')
self.skip_to_nonblank()
- tname='element'
- else: tname=last_tablename
+ tname = 'element'
+ else: tname = last_tablename
while tname != tablename:
self.skipto('@@@@@')
- tname=self.next_table_TOUGH2()
+ tname = self.next_table_TOUGH2()
def skip_to_table_TOUGHplus(self,tablename,last_tablename,nelt_tables):
if last_tablename is None:
self.skipto('=====',0)
self.skip_to_nonblank()
- tname='element'
- nelt_tables=0
- else: tname=last_tablename
+ tname = 'element'
+ nelt_tables = 0
+ else: tname = last_tablename
while tname != tablename:
- if tname=='primary': keyword='_____'
- else: keyword='@@@@@'
+ if tname == 'primary': keyword='_____'
+ else: keyword = '@@@@@'
self.skipto(keyword,0)
- tname=self.next_table_TOUGHplus()
- if tname=='element':
- nelt_tables+=1
- tname+=str(nelt_tables)
+ tname = self.next_table_TOUGHplus()
+ if tname == 'element':
+ nelt_tables += 1
+ tname += str(nelt_tables)
def start_of_values(self,line):
- """Returns start index of values in a table line. Characters before this start index are taken to contain
- the key(s) and row index."""
+ """Returns start index of values in a table line. Characters before
+ this start index are taken to contain the key(s) and row index."""
pt = line.find('.')
if pt >= 2:
nextpt = line.find('.', pt + 1)
@@ -586,28 +598,32 @@ def start_of_values(self,line):
if pos > 0 : return pos + 1
return None
- def key_positions(self,line,nkeys):
- """Returns detected positions of keys in the start of a table line. This works on the assumption
- that key values must have a digit present in their last character. It searches backwards from the
- end of the position of INDEX (or IND.) in the header line- sometimes keys can overlap into this area."""
- keylength=5
- keypos=[]
- pos=len(line)-1
- while line[pos]==' ': pos-=1
- while line[pos] != ' ': pos-=1
+ def key_positions(self, line, nkeys):
+ """Returns detected positions of keys in the start of a table line.
+ This works on the assumption that key values must have a digit
+ present in their last character. It searches backwards from
+ the end of the position of INDEX (or IND.) in the header line-
+ sometimes keys can overlap into this area. """
+ keylength = 5
+ keypos = []
+ pos = len(line)-1
+ while line[pos] == ' ': pos -= 1
+ while line[pos] != ' ': pos -= 1
for k in range(nkeys):
- while not line[pos].isdigit() and pos>=keylength: pos-=1
- pos-=(keylength-1)
- if valid_blockname(line[pos:pos+keylength]):
+ while not line[pos].isdigit() and pos >= keylength:
+ pos -= 1
+ pos -= (keylength - 1)
+ if valid_blockname(line[pos: pos + keylength]):
keypos.append(pos)
- pos-=1
+ pos -= 1
else: return None
keypos.reverse()
if len(keypos) != nkeys: return None
else: return keypos
def parse_table_header_AUTOUGH2(self):
- """Parses table header line for AUTOUGH2, returning the number of keys and the column names."""
+ """Parses table header line for AUTOUGH2, returning the number of keys
+ and the column names."""
cols = []
headline = self.readline()
headstrs = headline.strip().split()
@@ -630,16 +646,16 @@ def setup_table_AUTOUGH2(self,tablename):
# Double-check number of columns:
nvalues = len([s for s in line[start:].strip().split()])
if (len(cols) == nvalues):
- keypos = self.key_positions(line[:start],nkeys)
+ keypos = self.key_positions(line[:start], nkeys)
if keypos:
# determine row names:
while line[1:6] != keyword:
- keyval = [fix_blockname(line[kp:kp+5]) for kp in keypos]
- if len(keyval)>1: keyval = tuple(keyval)
+ keyval = [fix_blockname(line[kp: kp + 5]) for kp in keypos]
+ if len(keyval) > 1: keyval = tuple(keyval)
else: keyval = keyval[0]
rows.append(keyval)
line = self.readline()
- row_format = {'key':keypos, 'values':[start]}
+ row_format = {'key': keypos, 'values': [start]}
allow_rev = tablename == 'connection'
self._table[tablename] = listingtable(cols, rows, row_format, num_keys = nkeys,
allow_reverse_keys = allow_rev)
@@ -738,8 +754,10 @@ def is_separator(line): return len(line)>lsep and line[1:lsep+1] == line[1]*lsep
else: keyval = keyval[0]
indexstr = line[index_pos[0]:index_pos[1]]
try: index = int(indexstr) - 1
- except ValueError: index += 1 # to handle overflow (****) in index field: assume indices continue
- rowdict[index] = (count,keyval) # use a dictionary to deal with duplicate row indices (TOUGH2_MP)
+ # To handle overflow (****) in index field: assume indices continue:
+ except ValueError: index += 1
+ # Use a dictionary to deal with duplicate row indices (TOUGH2_MP):
+ rowdict[index] = (count,keyval)
if len(line.strip()) > len(longest_line): longest_line = line
pos = self._file.tell()
last_count = count
@@ -766,14 +784,15 @@ def is_separator(line): return len(line)>lsep and line[1:lsep+1] == line[1]*lsep
for i in range(internal_header_skiplines): line,count = count_read(count)
skiplines.append(count - last_count - 1)
indices = sorted(rowdict.keys())
- row_line=[rowdict[index][0] for index in indices]
- rows=[rowdict[index][1] for index in indices]
- numpos=self.parse_table_line(longest_line,start)
- row_format={'key':keypos,'index':keypos[-1]+5,'values':numpos}
- allow_rev=tablename=='connection'
+ row_line = [rowdict[index][0] for index in indices]
+ rows = [rowdict[index][1] for index in indices]
+ numpos = self.parse_table_line(longest_line,start)
+ row_format = {'key': keypos, 'index': keypos[-1] + 5, 'values': numpos}
+ allow_rev = tablename == 'connection'
self._table[tablename] = listingtable(cols, rows, row_format, row_line, num_keys = nkeys,
- allow_reverse_keys = allow_rev,
- header_skiplines = header_skiplines, skiplines = skiplines)
+ allow_reverse_keys = allow_rev,
+ header_skiplines = header_skiplines,
+ skiplines = skiplines)
self._tablenames.append(tablename)
else: raise Exception('Error parsing '+tablename+' table keys: table not created.')
@@ -831,12 +850,12 @@ def next_tablename(self, tablename):
else: return None
def read_tables_AUTOUGH2(self):
- tablename='element'
+ tablename = 'element'
while tablename:
self.read_header()
if tablename in self.skip_tables: self.skip_table(tablename)
else: self.read_table(tablename)
- tablename=self.next_table()
+ tablename = self.next_table()
def read_tables_TOUGH2(self):
tablename = 'element'
@@ -847,21 +866,22 @@ def read_tables_TOUGH2(self):
elif tablename in self._table: self.read_table(tablename)
else: # tables not present at first time step
next_tablename = self.next_tablename(last_tablename)
- if next_tablename: self.skip_to_table(next_tablename, last_tablename, 1)
+ if next_tablename:
+ self.skip_to_table(next_tablename, last_tablename, 1)
last_tablename = tablename
tablename = self.next_table()
def read_tables_TOUGHplus(self):
tablename='element'
self.read_header() # only one header at each time
- nelt_tables=0
+ nelt_tables = 0
while tablename:
if tablename in self.skip_tables: self.skip_table(tablename)
else: self.read_table(tablename)
- tablename=self.next_table()
- if tablename=='element':
- nelt_tables+=1
- tablename+=str(nelt_tables)
+ tablename = self.next_table()
+ if tablename == 'element':
+ nelt_tables += 1
+ tablename += str(nelt_tables)
def read_table_AUTOUGH2(self,tablename):
fmt = self._table[tablename].row_format
@@ -886,14 +906,15 @@ def skip_table_AUTOUGH2(self,tablename):
self._file.readline()
def read_table_line_AUTOUGH2(self,line,num_columns=None,fmt=None):
- start=fmt['values'][0]
- vals=[fortran_float(s) for s in line[start:].strip().split()]
+ start = fmt['values'][0]
+ vals = [fortran_float(s) for s in line[start:].strip().split()]
return vals
def read_table_line_TOUGH2(self,line,num_columns,fmt):
"""Reads values from a line in a TOUGH2 listing, given the number of columns, and format."""
- nvals=len(fmt['values'])-1
- return [fortran_float(line[fmt['values'][i]:fmt['values'][i+1]]) for i in range(nvals)]+[0.0]*(num_columns-nvals)
+ nvals = len(fmt['values']) - 1
+ return [fortran_float(line[fmt['values'][i]: fmt['values'][i+1]])
+ for i in range(nvals)] + [0.0] * (num_columns - nvals)
def read_table_TOUGH2(self,tablename):
table = self._table[tablename]
@@ -1050,8 +1071,9 @@ def datetime_array(t): return np.array([start_datetime + timedelta(0, s) for s i
return result
def get_reductions(self):
- """Returns a list of time step indices at which the time step is reduced, and the blocks at which the maximum
- residual occurred prior to the reduction."""
+ """Returns a list of time step indices at which the time step is
+ reduced, and the blocks at which the maximum residual occurred
+ prior to the reduction."""
self.rewind()
line, lastline = '', ''
keyword = "+++++++++ REDUCE TIME STEP"
@@ -1080,31 +1102,39 @@ def get_reductions(self):
return rl
reductions=property(get_reductions)
- def get_difference(self,indexa=None,indexb=None):
- """Returns dictionary of maximum differences, and locations of difference, of all element table properties between
- two sets of results. If both indexa and indexb are provided, the result is the difference between these two result indices.
- If only one index is given, the result is the difference between the given index and the one before that.
- If neither are given, the result is the difference between the last and penultimate sets of results."""
+ def get_difference(self, indexa = None, indexb = None):
+ """Returns dictionary of maximum differences, and locations of
+ difference, of all element table properties between two sets
+ of results. If both indexa and indexb are provided, the
+ result is the difference between these two result indices. If
+ only one index is given, the result is the difference between
+ the given index and the one before that. If neither are
+ given, the result is the difference between the last and
+ penultimate sets of results.
+ """
from copy import deepcopy
- tablename='element'
+ tablename = 'element'
if indexa is None: self.last()
else: self.set_index(indexa)
- results2=deepcopy(self._table[tablename])
+ results2 = deepcopy(self._table[tablename])
if indexb is None: self.prev()
else: self.set_index(indexb)
- results1=self._table[tablename]
- cvg={}
+ results1 = self._table[tablename]
+ cvg = {}
for name in results1.column_name:
- iblk=np.argmax(abs(results2[name]-results1[name]))
- blkname=results1.row_name[iblk]
- diff=results2[name][iblk]-results1[name][iblk]
- cvg[name]=(diff,blkname)
+ iblk = np.argmax(abs(results2[name] - results1[name]))
+ blkname = results1.row_name[iblk]
+ diff = results2[name][iblk] - results1[name][iblk]
+ cvg[name] = (diff,blkname)
return cvg
- convergence=property(get_difference)
+ convergence = property(get_difference)
- def get_vtk_data(self, geo, grid = None, flows = False, flux_matrix = None, geo_matches = True, blockmap = {}):
- """Returns dictionary of VTK data arrays from listing file at current time. If flows is True, average flux vectors
- are also calculated from connection data at the block centres."""
+ def get_vtk_data(self, geo, grid = None, flows = False, flux_matrix = None,
+ geo_matches = True, blockmap = {}):
+ """Returns dictionary of VTK data arrays from listing file at current
+ time. If flows is True, average flux vectors are also
+ calculated from connection data at the block centres.
+ """
from vtk import vtkFloatArray
natm = geo.num_atmosphere_blocks
nele = geo.num_underground_blocks
@@ -1115,7 +1145,8 @@ def get_vtk_data(self, geo, grid = None, flows = False, flux_matrix = None, geo_
flownames = []
def is_flowname(name):
name = name.lower()
- return name.startswith('flo') or name.endswith('flo') or name.endswith('flow') or name.endswith('veloc')
+ return name.startswith('flo') or name.endswith('flo') or \
+ name.endswith('flow') or name.endswith('veloc')
if flows:
if flux_matrix is None: flux_matrix = grid.flux_matrix(geo, blockmap)
flownames = [name for name in self.connection.column_name if is_flowname(name)]
@@ -1135,7 +1166,8 @@ def mname(blk): return blockmap[blk] if blk in blockmap else blk
array.SetNumberOfComponents(1)
array.SetNumberOfValues(array_length[array_type])
for tablename in elt_tablenames:
- if geo_matches: array_data[array_type][name] = self._table[tablename][name][natm:] # faster
+ if geo_matches:
+ array_data[array_type][name] = self._table[tablename][name][natm:] # faster
else: # more flexible
array_data[array_type][name] = \
np.array([self._table[tablename][mname(blk)][name]
@@ -1150,13 +1182,19 @@ def mname(blk): return blockmap[blk] if blk in blockmap else blk
arrays[array_type][name].SetValue(iblk, data[iblk])
return arrays
- def write_vtk(self, geo, filename, grid = None, indices = None, flows = False, wells = False, start_time = 0.0,
- time_unit = 's', flux_matrix = None, blockmap = {}, surface_snap = 0.1):
- """Writes VTK files for a vtkUnstructuredGrid object corresponding to the grid in 3D with the listing data,
- with the specified filename, for visualisation with VTK. A t2grid can optionally be specified, to include rock type
- data as well. A list of the required time indices can optionally be specified. If a grid is specified, flows is True,
- and connection data are present in the listing file, approximate average flux vectors are also calculated at the
- block centres from the connection data."""
+ def write_vtk(self, geo, filename, grid = None, indices = None, flows = False,
+ wells = False, start_time = 0.0, time_unit = 's',
+ flux_matrix = None, blockmap = {}, surface_snap = 0.1):
+ """Writes VTK files for a vtkUnstructuredGrid object corresponding to
+ the grid in 3D with the listing data, with the specified
+ filename, for visualisation with VTK. A t2grid can optionally
+ be specified, to include rock type data as well. A list of
+ the required time indices can optionally be specified. If a
+ grid is specified, flows is True, and connection data are
+ present in the listing file, approximate average flux vectors
+ are also calculated at the block centres from the connection
+ data.
+ """
from vtk import vtkXMLUnstructuredGridWriter
from os.path import splitext
base, ext = splitext(filename)
@@ -1165,7 +1203,8 @@ def write_vtk(self, geo, filename, grid = None, indices = None, flows = False, w
doflows = False
if flows and (self.connection is not None):
if grid is None:
- raise Exception("t2listing.write_vtk(): if flows == True, a t2grid object must be specified.")
+ raise Exception("t2listing.write_vtk(): if flows == True, " +
+ " a t2grid object must be specified.")
else:
if geo_matches or len(blockmap) > 0: doflows = True
else:
@@ -1186,7 +1225,7 @@ def write_vtk(self, geo, filename, grid = None, indices = None, flows = False, w
collection = pvd.createElement('Collection')
initial_index = self.index
if indices is None: indices = range(self.num_fulltimes)
- timescales = {'s':1.0,'h':3600.,'d':3600.*24,'y':3600.*24*365.25}
+ timescales = {'s': 1.0, 'h': 3600., 'd': 3600.*24, 'y': 3600. * 24 * 365.25}
if time_unit in timescales: timescale = timescales[time_unit]
else: timescale = 1.0
writer = vtkXMLUnstructuredGridWriter()
@@ -1213,7 +1252,7 @@ def write_vtk(self, geo, filename, grid = None, indices = None, flows = False, w
pvdfile.close()
self.index = initial_index
- def add_side_recharge(self,geo,dat):
+ def add_side_recharge(self, geo, dat):
"""Adds side recharge generators to a TOUGH2 data object for a production run,
calculated according to the final results in the listing. These generators represent side
inflows due to pressure changes in the blocks on the model's horizontal boundaries.
@@ -1222,104 +1261,116 @@ def add_side_recharge(self,geo,dat):
from IAPWS97 import cowat,visc
from geometry import line_projection
from t2data import t2generator
- initial_index=self.index
- keyword={'AUTOUGH2':{'P':'Pressure','T':'Temperature'},'TOUGH2':{'P':'P','T':'T'},
- 'TOUGH2_MP':{'P':'P','T':'T'},'TOUGH+':{'P':'Pressure','T':'Temperature'}}
+ initial_index = self.index
+ keyword = {'AUTOUGH2':{'P': 'Pressure', 'T': 'Temperature'},
+ 'TOUGH2': {'P': 'P', 'T': 'T'},
+ 'TOUGH2_MP': {'P': 'P', 'T': 'T'},
+ 'TOUGH+' : {'P': 'Pressure', 'T': 'Temperature'}}
self.last()
- bdy_nodes=geo.boundary_nodes
+ bdy_nodes = geo.boundary_nodes
for blk in dat.grid.blocklist[geo.num_atmosphere_blocks:]:
- colname=geo.column_name(blk.name)
+ colname = geo.column_name(blk.name)
if colname in geo.column:
- col=geo.column[colname]
- if col.num_neighbours0: nkeystr=str(self.num_keys)
- else: nkeystr='summed'
- return 'History results for '+nkeystr+' '+self.keytype+'s at '+str(self.num_times)+' times'
-
- def __getitem__(self,key):
- """Returns results for a given key value, in the form of a dictionary with keys corresponding to the
- column names. Each value is an array of history results for that key and column name- unless the key also
- has the time appended as its third element, in which case the dictionary values are just single floating
- point values for each column."""
- if self.num_rows>0:
- if self._nkeys>0:
- if not isinstance(key,tuple): key=(key,)
+ if self._nkeys > 0: nkeystr = str(self.num_keys)
+ else: nkeystr = 'summed'
+ return 'History results for ' + nkeystr + ' ' + self.keytype + \
+ 's at ' + str(self.num_times) + ' times'
+
+ def __getitem__(self, key):
+ """Returns results for a given key value, in the form of a dictionary
+ with keys corresponding to the column names. Each value is an
+ array of history results for that key and column name- unless
+ the key also has the time appended as its third element, in
+ which case the dictionary values are just single floating
+ point values for each column.
+ """
+ if self.num_rows > 0:
+ if self._nkeys > 0:
+ if not isinstance(key, tuple): key = (key,)
if key in self.keys:
- keydata=self._data[self._keyrows[key]]
- return dict([(colname,keydata[:,icol+1]) for icol,colname in enumerate(self.column_name)])
+ keydata = self._data[self._keyrows[key]]
+ return dict([(colname, keydata[:,icol+1]) for
+ icol, colname in enumerate(self.column_name)])
elif key in self._row:
- row=self._data[self._row[key]]
- return dict([(colname,row[icol+1]) for icol,colname in enumerate(self.column_name)])
+ row = self._data[self._row[key]]
+ return dict([(colname,row[icol+1]) for
+ icol, colname in enumerate(self.column_name)])
else: return None
else: # no keys (e.g. TOUGH+ COFT/GOFT)
try:
- icol=self.column_name.index(key)
- return self._data[:,icol+1]
+ icol = self.column_name.index(key)
+ return self._data[:, icol + 1]
except ValueError: return None
else: return None
- def read(self,filename):
+ def read(self, filename):
"""Reads contents of file(s) and stores in appropriate data structures."""
- self._rowindex=0
+ self._rowindex = 0
from glob import glob
- files=glob(filename)
- configured=False
- for i,fname in enumerate(files):
- self._file=open(fname,'rU')
- header=self._file.readline()
+ files = glob(filename)
+ configured = False
+ for i, fname in enumerate(files):
+ self._file = open(fname,'rU')
+ header = self._file.readline()
if header:
if not configured:
self.detect_simulator(header)
@@ -1329,199 +1380,214 @@ def read(self,filename):
self._file.close()
self.finalize_data()
- def detect_simulator(self,header):
+ def detect_simulator(self, header):
"""Detects simulator (TOUGH2 or TOUGH2_MP) from header line."""
if 'OFT' in header: self.simulator='TOUGH2_MP'
elif header.startswith('Time ['): self.simulator='TOUGH+'
else: self.simulator='TOUGH2'
- internal_fns=['setup_headers','read_data']
- simname=self.simulator.replace('+','plus')
+ internal_fns = ['setup_headers','read_data']
+ simname = self.simulator.replace('+', 'plus')
for fname in internal_fns:
- fname_sim=fname+'_'+simname
- setattr(self,fname,getattr(self,fname_sim))
-
- def setup_headers_TOUGH2(self,filename,header):
- """Sets up keys and column headings from given filename and header line, for TOUGH2 output."""
- if filename.endswith('OFT') and len(filename)>=4: self.type=filename[-4:].strip()
- else: self.type=None
- self.time_index=1
- self._nkeys=1
- items=header.strip().split(',')
- if items[-1]=='': del items[-1] # often an extra comma on the end of the lines
- items=items[3:]
- int_index=0
- for i,item in enumerate(items):
+ fname_sim = fname + '_' + simname
+ setattr(self, fname, getattr(self, fname_sim))
+
+ def setup_headers_TOUGH2(self, filename, header):
+ """Sets up keys and column headings from given filename and header
+ line, for TOUGH2 output."""
+ if filename.endswith('OFT') and len(filename) >= 4:
+ self.type = filename[-4:].strip()
+ else: self.type = None
+ self.time_index = 1
+ self._nkeys = 1
+ items = header.strip().split(',')
+ if items[-1] == '': del items[-1] # often an extra comma on the end of the lines
+ items = items[3:]
+ int_index = 0
+ for i, item in enumerate(items):
try:
- int_item=int(item)
- int_index=i
+ int_item = int(item)
+ int_index = i
break
except: pass
- if int_index==0: ncols=len(items)
- else: ncols=int_index
+ if int_index == 0: ncols = len(items)
+ else: ncols = int_index
self.column_name=range(ncols)
- def setup_headers_TOUGH2_MP(self,filename,header):
- """Sets up keys and column headings from given filename and header line, for TOUGH2_MP output."""
- headers=header.strip().split()
- self.type=headers[0] # FOFT, COFT or GOFT
- time_header=[h for h in headers if h.lower().startswith('time')][0]
- self.time_index=headers.index(time_header)
- if self.type=='FOFT':
- self.key_index=self.time_index-1
- self._nkeys=1
+ def setup_headers_TOUGH2_MP(self, filename, header):
+ """Sets up keys and column headings from given filename and header
+ line, for TOUGH2_MP output."""
+ headers = header.strip().split()
+ self.type = headers[0] # FOFT, COFT or GOFT
+ time_header = [h for h in headers if h.lower().startswith('time')][0]
+ self.time_index = headers.index(time_header)
+ if self.type == 'FOFT':
+ self.key_index = self.time_index - 1
+ self._nkeys = 1
else: # COFT or GOFT
- self.key_index=self.time_index+1
- self._nkeys=2
- self.key_name=headers[self.key_index:self.key_index+self._nkeys]
- prepend_titles, append_titles=['GAS','GENERATION'],['flow']
- startcol=self._nkeys+2
- cols=[]
- i=startcol
- while i<=len(headers)-1:
- title=headers[i]
+ self.key_index = self.time_index + 1
+ self._nkeys = 2
+ self.key_name = headers[self.key_index: self.key_index+self._nkeys]
+ prepend_titles, append_titles = ['GAS','GENERATION'], ['flow']
+ startcol = self._nkeys + 2
+ cols = []
+ i = startcol
+ while i <= len(headers) - 1:
+ title = headers[i]
if title in prepend_titles:
- cols.append(title+' '+headers[i+1])
- i+=1
+ cols.append(title + ' ' + headers[i+1])
+ i += 1
elif title in append_titles:
- cols[-1]+=' '+title
+ cols[-1] += ' ' + title
else: cols.append(title)
- i+=1
- self.column_name=cols
- self.col_start=[header.index(colname) for colname in self.column_name]
- self.key_start=[header.index(key) for key in self.key_name]
- self.time_pos=[header.index(time_header)]
- if self.type=='FOFT':
+ i += 1
+ self.column_name = cols
+ self.col_start = [header.index(colname) for colname in self.column_name]
+ self.key_start = [header.index(key) for key in self.key_name]
+ self.time_pos = [header.index(time_header)]
+ if self.type == 'FOFT':
self.key_start.append(self.time_pos[0])
self.time_pos.append(self.col_start[0])
else:
self.key_start.append(self.col_start[0])
self.time_pos.append(self.key_start[0])
- def setup_headers_TOUGHplus(self,filename,header):
- """Sets up keys and column headings from given filename and header line, for TOUGH+ output."""
- headers=header.strip().split('-')
- if filename.endswith('OFT') and len(filename)>=4: self.type=filename[-4:].strip()
+ def setup_headers_TOUGHplus(self, filename, header):
+ """Sets up keys and column headings from given filename and header
+ line, for TOUGH+ output."""
+ headers = header.strip().split('-')
+ if filename.endswith('OFT') and len(filename) >= 4:
+ self.type = filename[-4:].strip()
elif '_Time_Series' in filename:
- filetype={'Elem_Time_Series':'FOFT','Conx_Time_Series':'COFT','SS_Time_Series':'GOFT'}
+ filetype = {'Elem_Time_Series': 'FOFT', 'Conx_Time_Series': 'COFT',
+ 'SS_Time_Series': 'GOFT'}
for key in filetype.keys():
if key in filename:
- self.type=filetype[key]
+ self.type = filetype[key]
break
- else: self.type=None
- if self.type=='FOFT':
- self.time_index=1
- self._nkeys=1
+ else: self.type = None
+ if self.type == 'FOFT':
+ self.time_index = 1
+ self._nkeys = 1
else:
- self.time_index=0
- self._nkeys=0
- cols=headers[self._nkeys+1:]
+ self.time_index = 0
+ self._nkeys = 0
+ cols = headers[self._nkeys + 1:]
from re import sub
- cols=[sub('\[.*\]','',col).strip() for col in cols] # remove units
- self.column_name=cols
+ cols = [sub('\[.*\]','',col).strip() for col in cols] # remove units
+ self.column_name = cols
- def read_data_TOUGH2(self,configured):
+ def read_data_TOUGH2(self, configured):
"""Reads in the data, for TOUGH2 output."""
self._file.seek(0)
- lines=self._file.readlines()
+ lines = self._file.readlines()
for line in lines:
- items=line.strip().split(',')
- if items[-1]=='': del items[-1]
+ items = line.strip().split(',')
+ if items[-1] == '': del items[-1]
time_index = fortran_int(items.pop(0))
- time=float(items.pop(0))
+ time = float(items.pop(0))
self.times.append(time)
- nc1=self.num_columns+1
- nsets=len(items)/nc1
+ nc1 = self.num_columns + 1
+ nsets = len(items) / nc1
for i in range(nsets):
- setvals=items[i*nc1:(i+1)*nc1]
+ setvals = items[i*nc1: (i+1)*nc1]
key = (fortran_int(setvals[0]),)
- vals=[fortran_float(val) for val in setvals[1:]]
+ vals = [fortran_float(val) for val in setvals[1:]]
self.row_name.append(key+(time,))
if not key in self.keys:
- self._keyrows[key]=[]
+ self._keyrows[key] = []
self.keys.append(key)
self._keyrows[key].append(self._rowindex)
self._data.append([time]+vals)
- self._rowindex+=1
+ self._rowindex += 1
- def read_data_TOUGH2_MP(self,configured):
+ def read_data_TOUGH2_MP(self, configured):
"""Reads in the data, for TOUGH2_MP output."""
def get_key(line):
- return tuple([fix_blockname(line[self.key_start[i]:self.key_start[i+1]].rstrip()) for i in range(self._nkeys)])
- def get_time(line): return fortran_float(line[self.time_pos[0]:self.time_pos[1]])
+ return tuple([fix_blockname(line[self.key_start[i]:
+ self.key_start[i+1]].rstrip()) for
+ i in range(self._nkeys)])
+ def get_time(line):
+ return fortran_float(line[self.time_pos[0]:self.time_pos[1]])
def get_vals(line):
- start=self.col_start+[len(line)] # allow for lines of different lengths
- return [fortran_float(line[start[i]:start[i+1]]) for i in range(self.num_columns)]
- lines=self._file.readlines()
- first_key=None
+ start = self.col_start + [len(line)] # allow for lines of different lengths
+ return [fortran_float(line[start[i]: start[i+1]]) for
+ i in range(self.num_columns)]
+ lines = self._file.readlines()
+ first_key = None
from copy import copy
- otherfile_keys=copy(self.keys)
+ otherfile_keys = copy(self.keys)
for line in lines:
if line.strip():
- time=get_time(line)
- key=get_key(line)
- if not first_key: first_key=key
- vals=get_vals(line)
- rowname=key+(time,)
+ time = get_time(line)
+ key = get_key(line)
+ if not first_key: first_key = key
+ vals = get_vals(line)
+ rowname = key + (time,)
if not (key in otherfile_keys):
if key in self._keyrows:
- keyrows=self._keyrows[key]
- keytimes=[self._data[irow][0] for irow in self._keyrows[key]]
- keyrownames=[self.row_name[irow] for irow in self._keyrows[key]]
- newtime=not (time in keytimes)
- newrowname=not (rowname in keyrownames)
- else: newtime,newrowname=True,True
- if (not configured) and (key==first_key) and newtime: self.times.append(time)
+ keyrows = self._keyrows[key]
+ keytimes = [self._data[irow][0] for irow in self._keyrows[key]]
+ keyrownames = [self.row_name[irow] for irow in self._keyrows[key]]
+ newtime = not (time in keytimes)
+ newrowname = not (rowname in keyrownames)
+ else: newtime, newrowname = True, True
+ if (not configured) and (key == first_key) and newtime:
+ self.times.append(time)
if newrowname:
self.row_name.append(rowname)
if not key in self.keys:
- self._keyrows[key]=[]
+ self._keyrows[key] = []
self.keys.append(key)
self._keyrows[key].append(self._rowindex)
- self._data.append([time]+vals)
- self._rowindex+=1
+ self._data.append([time] + vals)
+ self._rowindex += 1
- def read_data_TOUGHplus(self,configured):
+ def read_data_TOUGHplus(self, configured):
"""Reads in the data, for TOUGH+ output."""
- lines=self._file.readlines()
- if self.type=='FOFT':
+ lines = self._file.readlines()
+ if self.type == 'FOFT':
for line in lines:
- items=line.strip().split(',')
- time=fortran_float(items[1])
+ items = line.strip().split(',')
+ time = fortran_float(items[1])
self.times.append(time)
- nc1=self.num_columns+1
- nsets=(len(items)-2)/nc1
+ nc1 = self.num_columns + 1
+ nsets = (len(items)-2) / nc1
for i in range(nsets):
- setvals=items[2+i*nc1:2+(i+1)*nc1]
+ setvals = items[2+i*nc1: 2+(i+1)*nc1]
key = (fortran_int(setvals[0]),)
- if key[0]>0:
- vals=[fortran_float(val) for val in setvals[1:]]
- self.row_name.append(key+(time,))
+ if key[0] > 0:
+ vals = [fortran_float(val) for val in setvals[1:]]
+ self.row_name.append(key + (time,))
if not key in self.keys:
- self._keyrows[key]=[]
+ self._keyrows[key] = []
self.keys.append(key)
self._keyrows[key].append(self._rowindex)
- self._data.append([time]+vals)
- self._rowindex+=1
+ self._data.append([time] + vals)
+ self._rowindex += 1
else:
for line in lines:
- vals=[fortran_float(val) for val in line.strip().split()]
- time=vals.pop(0)
+ vals = [fortran_float(val) for val in line.strip().split()]
+ time = vals.pop(0)
self.times.append(time)
- self._data.append([time]+vals)
+ self._data.append([time] + vals)
def finalize_data(self):
- self._data=np.array(self._data,float64)
- self.times=np.array(self.times,float64)
- self._row=dict([(r,i) for i,r in enumerate(self.row_name)])
+ self._data = np.array(self._data, float64)
+ self.times = np.array(self.times, float64)
+ self._row = dict([(r,i) for i,r in enumerate(self.row_name)])
class toughreact_tecplot(object):
- """Class for TOUGHREACT Tecplot output files. These work similarly to t2listing objects, but
- have just a single element table at each time. It is possible to navigate through time by
- using the next() and prev() functions to step through, or using the first() and last() functions to go to
- the start or end, or to set the index, step (model time step number) or time properties directly.
- When reading a toughreact_tecplot object from file it is necessary also to specify the block names
- (as these are not stored in the Tecplot file)."""
+ """Class for TOUGHREACT Tecplot output files. These work similarly to
+ t2listing objects, but have just a single element table at each
+ time. It is possible to navigate through time by using the next()
+ and prev() functions to step through, or using the first() and
+ last() functions to go to the start or end, or to set the index,
+ step (model time step number) or time properties directly. When
+ reading a toughreact_tecplot object from file it is necessary also
+ to specify the block names (as these are not stored in the Tecplot
+ file).
+ """
def __init__(self, filename, blocks):
self.filename = filename
self._file = open(filename, 'rU')
@@ -1621,14 +1687,17 @@ def setup_table(self, blocks):
if isinstance(blocks, mg.mulgrid): blocks = blocks.block_name_list
elif isinstance(blocks, t2g.t2grid): blocks = [blk.name for blk in blocks.blocklist]
if len(blocks) != self._num_blocks:
- raise Exception("Specified block name list is the wrong length for TOUGHREACT Tecplot file "+ self.filename)
+ raise Exception("Specified block name list is the wrong length for " +
+ "TOUGHREACT Tecplot file "+ self.filename)
self._file.seek(0)
line = self.skipto('VARIABLES')
if line is not None:
eqpos = line.find('=')
cols = [col.strip() for col in line[eqpos+1:].strip().split(',')[:-1]]
self.element = listingtable(cols, blocks, num_keys = 1)
- else: raise Exception("Could not find variable definitions for TOUGHREACT Tecplot file " + self.filename)
+ else:
+ raise Exception("Could not find variable definitions " +
+ "for TOUGHREACT Tecplot file " + self.filename)
def read_table_line(self, line):
"""Parses given string and returns an array of float values."""
@@ -1641,7 +1710,8 @@ def read_table(self):
self.element[i] = self.read_table_line(line)
def history(self, selection):
- """Returns time histories for specified selection of block names (or index) and column names."""
+ """Returns time histories for specified selection of block names (or
+ index) and column names."""
def ordered_selection(selection):
osel = []
@@ -1675,7 +1745,7 @@ def ordered_selection(selection):
hist[sel_index].append(vals[valindex])
self._index = old_index
result = [(self.times,np.array(h)) for sel_index,h in enumerate(hist)]
- if len(result)==1: result = result[0]
+ if len(result) == 1: result = result[0]
return result
def get_vtk_data(self, geo, grid = None, geo_matches = True, blockmap = {}):
@@ -1700,13 +1770,18 @@ def mname(blk): return blockmap[blk] if blk in blockmap else blk
for array_type,data_dict in array_data.items():
for name,data in data_dict.items():
- for iblk in range(nele): arrays[array_type][name].SetValue(iblk, data[iblk])
+ for iblk in range(nele):
+ arrays[array_type][name].SetValue(iblk, data[iblk])
return arrays
- def write_vtk(self, geo, filename, grid = None, indices = None, start_time = 0.0, time_unit = 's', blockmap = {}, surface_snap = 0.1):
- """Writes VTK files for a vtkUnstructuredGrid object corresponding to the grid in 3D with the Tecplot data,
- with the specified filename, for visualisation with VTK. A t2grid can optionally be specified, to include rock type
- data as well. A list of the required time indices can optionally be specified."""
+ def write_vtk(self, geo, filename, grid = None, indices = None, start_time = 0.0,
+ time_unit = 's', blockmap = {}, surface_snap = 0.1):
+ """Writes VTK files for a vtkUnstructuredGrid object corresponding to
+ the grid in 3D with the Tecplot data, with the specified
+ filename, for visualisation with VTK. A t2grid can optionally
+ be specified, to include rock type data as well. A list of
+ the required time indices can optionally be specified.
+ """
from vtk import vtkXMLUnstructuredGridWriter
from os.path import splitext
base, ext = splitext(filename)
@@ -1725,15 +1800,17 @@ def write_vtk(self, geo, filename, grid = None, indices = None, start_time = 0.0
initial_index = self.index
if indices is None: indices = range(self.num_times)
yr = 3600.*24*365.25
- timescales = {'s':1.0,'h':3600.,'d':3600.*24,'y':yr}
- if time_unit in timescales: timescale = timescales[time_unit]/yr # assumes Tecplot times are in years
+ timescales = {'s': 1.0, 'h': 3600., 'd': 3600. * 24, 'y': yr}
+ if time_unit in timescales:
+ timescale = timescales[time_unit] / yr # assumes Tecplot times are in years
else: timescale = 1.0
writer = vtkXMLUnstructuredGridWriter()
for i in indices:
self.index = i
t = start_time + self.time/timescale
- filename_time = base+'_'+str(i)+'.vtu'
- results_arrays = self.get_vtk_data(geo, grid, geo_matches=geo_matches, blockmap = blockmap)
+ filename_time = base + '_' + str(i) + '.vtu'
+ results_arrays = self.get_vtk_data(geo, grid, geo_matches = geo_matches,
+ blockmap = blockmap)
for array_type,array_dict in arrays.items():
array_dict.update(results_arrays[array_type])
vtu = geo.get_vtk_grid(arrays, surface_snap)