Source code for BETRS
################################################################################
# BETR-Research
#
# A framework to create spatially resolved multimedia fate- and transport
# models for chemical contaminants.
#
# Copyright (C) 2010 Harald von Waldow <hvwaldow@chem.ethz.ch>
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
################################################################################
from numpy import *
import sys
import pdb
import cPickle
import os
## these imports reload modules in each run, for making development
## workflow smoother
## adding the subdirectory "PY" to the module search path. "PY" contains
## model-specific code.
pydir=os.path.abspath(os.path.join(os.path.dirname(__file__),'PY'))
if not pydir in sys.path:
sys.path.insert(0,pydir)
import readinput
import globalz
import zvalues
import tempcorrect
import volumes
import processes
import mkflowD
import mkbigD
import emissions
import helpers
import solver
import output
import modifyparams
reload(readinput)
reload(globalz)
reload(zvalues)
reload(tempcorrect)
reload(volumes)
reload(processes)
reload(mkflowD)
reload(mkbigD)
reload(emissions)
reload(helpers)
reload(solver)
reload(output)
reload(modifyparams)
from readinput import *
from globalz import *
from zvalues import *
from tempcorrect import *
from volumes import *
from processes import *
from mkflowD import *
from mkbigD import *
from emissions import *
from helpers import *
from solver import *
from output import *
from modifyparams import *
##########################
[docs]class Model():
"""
The class :class:`Model` is the user\'s interface to BETR-Research.
:param int chemical: The ID of the chemical to be modelled (according to :file:`Chemicals/{chemdb}`)
:param string run: String describing a particular model run. Used for the construction of output file paths.
:param string chemdb: Filename of the chemicals database.
:param string seasonalparfile: Filename of the parameter-file for seasonally changing values.
:param string constantparfile: Filename of parameter-file for time-constant parameters.
:param string compartmentfile: Filename of file describing the compartment.
:param string flowdir: Name of the directory containing the convective flow descriptions.
:param string processfile: Filename of the list of processes to be considered.
:param string controlfile: Filename of the file containing options for this model run.
"""
def __init__(self,chemical, run='default',
chemdb='chemicals_default.txt',
seasonalparfile='seasonal_parameters_default.txt',
constantparfile='const_parameters_default.txt',
compartmentfile='compartments_default.txt',
flowdir='default',
processfile='processes_default.txt',
controlfile='control_default.txt'):
self.run=run
self.chemical=chemical
## read chemical
fn=os.path.join('Chemicals', chemdb)
self.update_chemical(chemical, fn)
## read environnmental parameters
fn1=os.path.join('Environment', constantparfile)
fn2=os.path.join('Environment', seasonalparfile)
self.update_env_params(fn1,fn2)
## read compartment descriptions
fn=os.path.join('Environment', compartmentfile)
self.update_compartments(fn)
## read flows of transport media
directory=os.path.join('Flows', flowdir)
self.update_flows(directory)
## read processlist
fn=os.path.join('Processes', processfile)
self.update_processlist(fn)
## read controlfile
fn=os.path.join('Control',controlfile)
self.update_control(fn)
## update model-characteristic class attributes
## nocells, nots, nocomps, matdim
self.update_char_numbers()
## construct time-averaged parameter array and flow-matrices
## parmean and flowdictmean
#self.average_parameters()
#self.average_flows()
## Calculate volumes of compartments and sub-compartments [self.vdict]
self.update_volumes() #
## update self.killidx, the vector of indices in 0:matdim indicating
## zero-volume compartments
self.mkkillidx_fromV()
## Calculate chemical properties in all cells at all timesteps
self.update_chempardict() #self.chempardict
## Calculate z-values [self.zdict]
self.update_zvalues()
## calculate D-values for intra-cell processes [self.Dproc]
self.update_dvalues_process()
## calculate D-values for flow in transport media [self.Dflow]
self.update_dvalues_flows()
## make 1/ZV available as list of diagonal sparse
## csr matrices [self.ZVinvlist]
self.update_zvinvlist()
# Construct the system matrices of D-values,
# in sparse "coo" format, one for each timestep.
# Compress matrices (delete 0 volume compartments)
# and remember positions of deleted rows/columns
# in self.killidx
self.update_bigDlist()
#################### end init ##############################################
## functions that read input files and construct or update the
## class attributes
## chemdict, par, compdict, flowdict, proclist, controldict
[docs] def update_chemical(self,chemno,chemdb):
"""reads the chemical to be modelled from the database of chemicals and updates Model.chemdict.
:param chemno: ID of the chemical in chemdb.
:param chemdb: path to the database of chemicals."""
self.chemdict = readChemicals([chemno], chemdb)[chemno]
def update_env_params(self,cpfile,spfile):
self.par = constructEnvironment(cpfile,spfile)
self.par = modparams(self)
def update_compartments(self, compfile):
self.compdict = readCompartments(compfile)
def update_flows(self, flowdir, mean=False):
self.flowdict = readFlows(self.compdict.keys(),flowdir)
def update_processlist(self, processfile):
self.proclist = readProcesses(self.compdict, processfile)
def update_control(self, controlfile):
self.controldict = readControl(controlfile)
############################################################################
## update of model-characteristic class attributes
def update_char_numbers(self):
self.nocells = self.par.shape[0] # number of cells
self.nots = self.par.shape[1] # number of timesteps in a period
self.nocomp=len(self.compdict) # number of compartments
self.matdim=self.nocells*self.nocomp # system-matrix dimension
############################################################################
## functions that operate on input parameters to
## create or update dictionaries:
## vdict:{<compartment>:{<subcompartment>:array((nocells, notimesteps))}}
## chempardict:{<compartment>:array((nocells,notimesteps),
## dtype=[k_reac,Kaw,Koa,Kow])
## zdict:{<compartment>:{<subcompartment>:array((nocells, notimesteps))}}
def update_volumes(self, mean=False):
if not mean:
self.vdict = Volumes(self.par,self.compdict).V
else:
self.vdict = Volumes(self.parmean,self.compdict).V
def update_chempardict(self, mean=False):
if not mean:
self.chempardict=tempcorrect(self.par, self.compdict, self.chemdict)
else:
self.chempardict=tempcorrect(self.parmean, self.compdict,\
self.chemdict)
def update_zvalues(self, mean=False):
if not mean:
self.zdict = Zvalues(self.par, self.compdict, self.chempardict).Z
else:
self.zdict=Zvalues(self.parmean, self.compdict, self.chempardict).Z
###########################################################################
## Calculate D-values
def update_dvalues_process(self):
self.procobj=process(self)
self.Dproc=self.procobj.getD()
def update_dvalues_flows(self):
self.Dflow=mkflowD(self)
############################################################################
## Construction of system matrices
[docs] def update_zvinvlist(self):
'''Make list of 1/ZV diagonal matrices (matdim,matdim)
in sparse csr-format'''
print('Making 1/ZV - matrices\n')
self.ZVinvlist=mkZVinv(self)
def mkkillidx_fromV(self):
#Calculating system-matrix indices indicating
#zero-volume compartments; setting self.killidx
killidx=[]
for c in self.compdict.keys():
zerovols=where(sum(self.vdict[c]['bulk'], axis=1)==0)[0]
killidx.extend(tocell(zerovols+1,c,self))
self.killidx=sort(array(killidx))
def mkkillidx_fromD(self):
#Calculating 0-volume compatment indices from Dvalue-matrices.
#Used for consistency check.
zerovols=[]
for m in self.bigDlist:
zerovols.append(findemptycells(m))
if not all([x==zerovols[0] for x in zerovols]):
sys.exit("update_bigDlist: detected inconsistent zero-volumes \
between timesteps. Aborting !")
return(array(zerovols[0]))
def _DtoK(self):
#Create mass rate-constant matrices from D-value matrices
for m in arange(0,len(self.bigDlist)):
self.bigDlist[m]= self.bigDlist[m]*(self.ZVinvlist[m])
[docs] def update_bigDlist(self):
'''Make system matrices ready for solving'''
self.bigDlist=mkbigD(self)
if self.controldict['dumpDmatrices'] in ['1','y','Y','yes','Yes','YES']:
fn=os.path.join('Output', self.run,'matrixdump.cpk')
if not os.path.exists(os.path.dirname(fn)):
os.mkdir(os.path.dirname(fn))
else:
if os.path.exists(fn):
print('Attention: overwriting %s\n') % (fn)
of=open(fn,'w')
cPickle.dump(self.bigDlist,of,protocol=2)
of.close()
print('wrote D-value - matrices to %s\n') % (fn)
## consistency check for killidx
print('Consistency check of zero-volume compartments: '),
kidx=self.mkkillidx_fromD()
if not all(kidx == self.killidx):
sys.exit('Detected inconsistencies between\
0-volumes and 0 rows/columns in D-matrices')
print(' passed\n')
## compress D-value list
for m in arange(0,len(self.bigDlist)):
self.bigDlist[m]=compressmat(self.bigDlist[m], self.killidx)
## transform bigDlist to mass rate constants (divide by ZV)
## and save in csr - sparse format
self._DtoK()
############################################################################
## Functions reading input that are not called in __init__
[docs] def update_emissions(self,emfile):
''' Constructs emission object from emissions-input-file'''
fn=os.path.join('Emissions', emfile)
self.emission = Emission(fn)
[docs] def update_solver(self, stepfile='stepping_default.txt',
solvparamfile='solvparams_default.txt',
initfile='initial_default.txt'):
'''
Updates or initializes the solver.
:param string stepfile: File describing the mapping
between timesteps and actual time
:param string solvparamfile: File containing the parameters
for the solver
:param string initfile: File containing the initial
contamination of the environments.
'''
stepfile=os.path.join('Solver', stepfile)
solvparamfile=os.path.join('Solver', solvparamfile)
initfile=os.path.join('Solver', initfile)
self.solver=Solver(self,stepfile,solvparamfile,initfile)
## Functions that solve the system
[docs] def solve_ss(self):
'''
Averages the all mass rate matrices for each season (month) to
get a time invariant system. Solves for the steady-state.
'''
## average D-value matrices
self.Davg = mean(self.bigDlist)
try:
q=self.emission.get_emission(0,self,type='array')
except AttributeError:
print("update_emissions(<filename>) has to be called before\n"\
+"solving for steady-state. Did nothing!\n")
return()
ss_res =dot(-inv(self.Davg.toarray()), q)
self.ss_res=expandvec(ss_res, self.killidx)
[docs] def solve_dyn(self, no_periods):
'''
Solves the dynamic model for *no_periods* (years).
'''
self.dyn_res=self.solver.solve_ode(no_periods)
reslist=[]
for resvec in self.dyn_res.T:
reslist.append(expandvec(resvec, self.killidx))
self.dyn_res=array(reslist).T
############################################################################
## output routines
[docs] def output_ss(self, filename='ss', units=['mol'], netcdf=False):
"""
:param string filename: Prefix of the output file(s). The output files
will be written to Output/<run>/<filename>_x.
:param list units: A list of strings to indicate the unit(s) of the
output. Possible list elements are 'mol', 'mol_per_m3', 'kg',
'kg_per_m3' and 'Pa'.
:param bool netcdf: Indicates whether netCDF-files are written or not.
"""
## filename without extension (automatically generated from filetype)
fn=os.path.join('Output',self.run,filename)
write_output_ss(self,fn,units,netcdf)
[docs] def output_dyn(self,filename='dyn',units=['mol'], netcdf=False):
"""
:param string filename: Prefix of the output file(s). The output files
will be written to Output/<run>/<filename>_x.
:param list units: A list of strings to indicate the unit(s) of the
output. Possible list elements are 'mol', 'mol_per_m3', 'kg',
'kg_per_m3' and 'Pa'.
:param bool netcdf: Indicates whether netCDF-files are written or not.
"""
## filename without extension (automatically generated from filetype)
fn=os.path.join('Output',self.run,filename)
write_output_dyn(self,fn,units,netcdf)
############################################################################
## accessors
[docs] def get_avgMatrix(self):
"""
Returns the averages system matrix of mass rate constants used
for the steady-state solution.
"""
print("Constructing full steady-state-matrix ...")
try:
avgmatrix=expandmatrix(self.Davg, self.killidx).toarray()
except AttributeError:
print("solve_ss() has to be called before acessing the averaged"\
+" system matrix.")
return()
print("OK")
return(avgmatrix)
[docs] def get_ssEmission(self):
"""
Returns the full emission vector used for steady-state
calculation in Compressed Sparse Column format.
"""
try:
ssem=self.emission.get_emission(0,self,'array')
except AttributeError:
print("update_emissions(<filename>) has to be called before"\
+" acessing the emission vector. Did nothing!")
return()
print("Expanding emission vector ...")
ssem=expandvec(ssem, self.killidx)
print("OK")
return(ssem)
## construct or update time-averages of parameters and flow-matrices
def average_parameters(self):
self.parmean=empty((self.par.shape[0]), dtype=self.par.dtype)
for field in self.par.dtype.names:
self.parmean[field]=mean(self.par[field], axis=1)
def average_flows(self):
self.flowdictmean={}
for item in self.flowdict.items():
meanmat=c_[item[1][:,0:2], mean(item[1][:,2:], axis=1)]
self.flowdictmean[item[0]]=meanmat
############################# END OF BETRS.py ##################################
