Skip to content

Commit

Permalink
fix(Mflist): allow None as a list entry (#765)
Browse files Browse the repository at this point in the history 
* fix(Mflist): allow None, 0, and -1 as a list entry

Mflist did not allow a stress_period_dictionary to have None, 0, or -1 as a valid entry for a stress period.  But MODFLOW and MT3D do allow this.  Mflist was refactored so that a stress period dictionary could be spd = {0: None, 1: -1, 1: 0} or other combinations on this.  This is necessary for the mt3d ssm package, but it also was necessary for loading a valid well file that might have ITMP = -1 for all stress periods.  This would be a strange case for MODFLOW, but ITMP values of -1 for the first stress period of a wel package is valid for MODFLOW.

* slight modification for nss==0 in first period

closes #743
  • Loading branch information
@langevin-usgs
langevin-usgs authored Dec 16, 2019
1 parent 8d35383 commit 51d8064
Show file tree
Hide file tree
Showing 5 changed files with 142 additions and 58 deletions.
71 changes: 45 additions & 26 deletions autotest/t060_test_lkt.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@
import flopy

# make the working directory
tpth = os.path.join('temp', 't057')
tpth = os.path.join('temp', 't060')
if not os.path.isdir(tpth):
os.makedirs(tpth)

Expand All @@ -25,7 +25,8 @@ def test_lkt_with_multispecies():
mtexe = 'mt3dusgs'

# Instantiate MODFLOW object in flopy
mf = flopy.modflow.Modflow(modelname=modelname, exe_name=mfexe, model_ws=modelpth, version='mfnwt')
mf = flopy.modflow.Modflow(modelname=modelname, exe_name=mfexe,
model_ws=modelpth, version='mfnwt')

Lx = 27500.0
Ly = 22000.0
Expand Down Expand Up @@ -58,8 +59,10 @@ def test_lkt_with_multispecies():
iprnwt = 1
ibotav = 1

nwt = flopy.modflow.ModflowNwt(mf, headtol=headtol, fluxtol=fluxtol, maxiterout=maxiterout,
thickfact=thickfact, linmeth=linmeth, iprnwt=iprnwt, ibotav=ibotav,
nwt = flopy.modflow.ModflowNwt(mf, headtol=headtol,
fluxtol=fluxtol, maxiterout=maxiterout,
thickfact=thickfact, linmeth=linmeth,
iprnwt=iprnwt, ibotav=ibotav,
options='SIMPLE')

## Instantiate discretization (DIS) package for MODFLOW-NWT
Expand Down Expand Up @@ -124,24 +127,25 @@ def test_lkt_with_multispecies():

# Create the discretization object
# itmuni = 4 (days); lenuni = 2 (meters)
dis = flopy.modflow.ModflowDis(mf, nlay, nrow, ncol, nper=1, delr=delr, delc=delc,
top=top1, botm=botm, laycbd=0, itmuni=4, lenuni=2,
steady=Steady, nstp=nstp, tsmult=tsmult, perlen=perlen)
dis = flopy.modflow.ModflowDis(mf, nlay, nrow, ncol, nper=1,
delr=delr, delc=delc,
top=top1, botm=botm, laycbd=0,
itmuni=4, lenuni=2,
steady=Steady, nstp=nstp, tsmult=tsmult,
perlen=perlen)

## Instantiate upstream weighting (UPW) flow package for MODFLOW-NWT

# UPW parameters
# UPW must be instantiated after DIS. Otherwise, during the mf.write_input() procedures,
# flopy will crash.


# First line of UPW input is: IUPWCB HDRY NPUPW IPHDRY
hdry = -1.e+30
iphdry = 0

# Next variables are: LAYTYP, LAYAVG, CHANI, LAYVKA, LAYWET
laytyp = [1, 1, 1] # >0: convertible
layavg = 0 # 0: harmonic mean
chani = 1.0 # >0: CHANI is the horizontal anisotropy for the entire layer
chani = 1.0 # >0: CHANI is the horizontal anis for entire layer
layvka = 0 # =0: indicates VKA is vertical hydraulic conductivity
laywet = 0 # Always set equal to zero in UPW package
hk = 3.172E-03
Expand All @@ -150,8 +154,10 @@ def test_lkt_with_multispecies():
ss = 0.00001
sy = 0.30

upw = flopy.modflow.ModflowUpw(mf, laytyp=laytyp, layavg=layavg, chani=chani, layvka=layvka,
laywet=laywet, ipakcb=53, hdry=hdry, iphdry=iphdry, hk=hk,
upw = flopy.modflow.ModflowUpw(mf, laytyp=laytyp, layavg=layavg,
chani=chani, layvka=layvka,
laywet=laywet, ipakcb=53, hdry=hdry,
iphdry=iphdry, hk=hk,
vka=vka, ss=ss, sy=sy)

## Instantiate basic (BAS or BA6) package for MODFLOW-NWT
Expand Down Expand Up @@ -482,21 +488,30 @@ def test_lkt_with_multispecies():
[0.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 0.0]]}

lak = flopy.modflow.ModflowLak(mf, nlakes=nlakes, ipakcb=ipakcb, theta=theta,
nssitr=nssitr, sscncr=sscncr, surfdep=surfdep,
stages=stages, stage_range=stage_range, lakarr=lkarr,
bdlknc=bdlknc, flux_data=flux_data, unit_number=16)
lak = flopy.modflow.ModflowLak(mf, nlakes=nlakes, ipakcb=ipakcb,
theta=theta,
nssitr=nssitr, sscncr=sscncr,
surfdep=surfdep,
stages=stages, stage_range=stage_range,
lakarr=lkarr,
bdlknc=bdlknc, flux_data=flux_data,
unit_number=16)

## Instantiate linkage with mass transport routing (LMT) package for MODFLOW-NWT (generates linker file)
## Instantiate linkage with mass transport routing (LMT) package
# for MODFLOW-NWT (generates linker file)

lmt = flopy.modflow.ModflowLmt(mf, output_file_name='lkttest.ftl', output_file_header='extended',
output_file_format='formatted', package_flows = ['lak'])
lmt = flopy.modflow.ModflowLmt(mf, output_file_name='lkttest.ftl',
output_file_header='extended',
output_file_format='formatted',
package_flows = ['lak'])


## Now work on MT3D-USGS file creation

mt = flopy.mt3d.Mt3dms(modflowmodel=mf, modelname=modelname, model_ws=modelpth,
version='mt3d-usgs', namefile_ext='mtnam', exe_name=mtexe,
mt = flopy.mt3d.Mt3dms(modflowmodel=mf, modelname=modelname,
model_ws=modelpth,
version='mt3d-usgs', namefile_ext='mtnam',
exe_name=mtexe,
ftlfilename='lkttest.ftl', ftlfree=True)

## Instantiate basic transport (BTN) package for MT3D-USGS
Expand All @@ -523,8 +538,10 @@ def test_lkt_with_multispecies():

btn = flopy.mt3d.Mt3dBtn(mt, lunit=lunit, ncomp=ncomp, mcomp=mcomp,
sconc=sconc, prsity=prsity, cinact=cinact,
laycon=laycon, thkmin=thkmin, nprs=nprs, nprobs=nprobs,
chkmas=True,nprmas=nprmas, perlen=perlen, dt0=dt0,
laycon=laycon, thkmin=thkmin, nprs=nprs,
nprobs=nprobs,
chkmas=True,nprmas=nprmas, perlen=perlen,
dt0=dt0,
nstp=nstp, tsmult=tsmult, mxstrn=mxstrn,
ttsmult=ttsmult, ttsmax=ttsmax, sconc2=sconc2)

Expand All @@ -538,7 +555,8 @@ def test_lkt_with_multispecies():
adv = flopy.mt3d.Mt3dAdv(mt, mixelm=mixelm, percel=percel, mxpart=mxpart,
nadvfd=nadvfd)

## Instantiate generalized conjugate gradient solver (GCG) package for MT3D-USGS
## Instantiate generalized conjugate gradient solver (GCG)
# package for MT3D-USGS

mxiter = 1
iter1 = 50
Expand Down Expand Up @@ -572,7 +590,8 @@ def test_lkt_with_multispecies():

lkt = flopy.mt3d.Mt3dLkt(mt, nlkinit=nlkinit, mxlkbc=mxlkbc, icbclk=icbclk,
ietlak=ietlak, coldlak=coldlak,
lk_stress_period_data=lkt_flux_data, coldlak2=coldlak2)
lk_stress_period_data=lkt_flux_data,
coldlak2=coldlak2)

mf.write_input()
mt.write_input()
Expand Down
103 changes: 79 additions & 24 deletions autotest/t068_test_ssm.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,37 +3,45 @@
"""

import os
import warnings
import flopy
import numpy as np

mf_exe_name = 'mf2005'
mt_exe_name = 'mt3dms'
v1 = flopy.which(mf_exe_name)
v2 = flopy.which(mt_exe_name)
run = True
if v1 is None or v2 is None:
run = False


def test_mt3d_ssm_with_nodata_in_1st_sp():
model_ws = os.path.join('.', 'temp', 't068')

nlay, nrow, ncol = 10, 10, 10

nlay, nrow, ncol = 3, 5, 5
perlen = np.zeros((10), dtype=np.float) + 10
nper = len(perlen)

ibound = np.ones((nlay,nrow,ncol), dtype=np.int)

botm = np.arange(-1,-11,-1)
botm = np.arange(-1,-4,-1)
top = 0.

# creating MODFLOW model

model_ws = 'data'
modelname = 'ssm_ex2'
model_ws = os.path.join('.', 'temp', 't068a')
modelname = 'model_mf'

mf = flopy.modflow.Modflow(modelname, model_ws=model_ws, version='mfnwt')
mf = flopy.modflow.Modflow(modelname, model_ws=model_ws,
exe_name=mf_exe_name)
dis = flopy.modflow.ModflowDis(mf, nlay=nlay, nrow=nrow, ncol=ncol,
perlen=perlen, nper=nper, botm=botm, top=top,
steady=False)

bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=top)
upw = flopy.modflow.ModflowUpw(mf, hk=100, vka=100, ss=0.00001, sy=0.1)
lpf = flopy.modflow.ModflowLpf(mf, hk=100, vka=100, ss=0.00001, sy=0.1)
oc = flopy.modflow.ModflowOc(mf)
nwt = flopy.modflow.ModflowNwt(mf)
pcg = flopy.modflow.ModflowPcg(mf)
lmt = flopy.modflow.ModflowLmt(mf)

# recharge
rchrate = {}
Expand All @@ -51,18 +59,18 @@ def test_mt3d_ssm_with_nodata_in_1st_sp():
# after the first stress period (this was crashing flopy
# version 3.2.13
ghb_data = {}
ghb_data[2] = [(4, 4, 4, 0.1, 1.5)]
ssm_data[2] = [(4, 4, 4, 1.0, itype['GHB'], 1.0, 100.0)]
ghb_data[5] = [(4, 4, 4, 0.25, 1.5)]
ssm_data[5] = [(4, 4, 4, 0.5, itype['GHB'], 0.5, 200.0)]
ssm_data[2] = [(nlay - 1, 4, 4, 1.0, itype['GHB'], 1.0, 100.0)]
ghb_data[2] = [(nlay - 1, 4, 4, 0.1, 1.5)]
ghb_data[5] = [(nlay - 1, 4, 4, 0.25, 1.5)]
ssm_data[5] = [(nlay - 1, 4, 4, 0.5, itype['GHB'], 0.5, 200.0)]

for k in range(nlay):
for i in range(nrow):
ghb_data[2].append((k, i, 0, 0.0, 100.0))
ssm_data[2].append((k, i, 0, 0.0, itype['GHB'], 0.0, 0.0))

ghb_data[5] = [(4, 4, 4, 0.25, 1.5)]
ssm_data[5] = [(4, 4, 4, 0.5, itype['GHB'], 0.5, 200.0)]
ghb_data[5] = [(nlay - 1, 4, 4, 0.25, 1.5)]
ssm_data[5] = [(nlay - 1, 4, 4, 0.5, itype['GHB'], 0.5, 200.0)]
for k in range(nlay):
for i in range(nrow):
ghb_data[5].append((k, i, 0, -0.5, 100.0))
Expand All @@ -72,7 +80,9 @@ def test_mt3d_ssm_with_nodata_in_1st_sp():
ghb = flopy.modflow.ModflowGhb(mf, stress_period_data=ghb_data)

# create MT3D-USGS model
mt = flopy.mt3d.Mt3dms(modflowmodel=mf, modelname=modelname, model_ws=model_ws)
modelname = 'model_mt'
mt = flopy.mt3d.Mt3dms(modflowmodel=mf, modelname=modelname,
model_ws=model_ws, exe_name='mt3dms')
btn = flopy.mt3d.Mt3dBtn(mt, sconc=0, ncomp=2, sconc2=50.0)
adv = flopy.mt3d.Mt3dAdv(mt)
ssm = flopy.mt3d.Mt3dSsm(mt, stress_period_data=ssm_data)
Expand All @@ -81,17 +91,62 @@ def test_mt3d_ssm_with_nodata_in_1st_sp():
# Write the output
mf.write_input()
mt.write_input()

# confirm that MT3D files exist
assert os.path.isfile(os.path.join(model_ws, '{}.{}'.format(mt.name, btn.extension[0]))) is True
assert os.path.isfile(os.path.join(model_ws, '{}.{}'.format(mt.name, adv.extension[0]))) is True
assert os.path.isfile(os.path.join(model_ws, '{}.{}'.format(mt.name, ssm.extension[0]))) is True
assert os.path.isfile(os.path.join(model_ws, '{}.{}'.format(mt.name, gcg.extension[0]))) is True


# run the models
if run:
success, buff = mf.run_model(report=True)
assert success, 'MODFLOW did not run'
success, buff = mt.run_model(report=True, normal_msg='Program completed.')
assert success, 'MT3D did not run'

model_ws2 = os.path.join('.', 'temp', 't068b')
mf2 = flopy.modflow.Modflow.load('model_mf.nam', model_ws=model_ws,
exe_name='mf2005')
mf2.change_model_ws(model_ws2)
mt2 = flopy.mt3d.Mt3dms.load('model_mt.nam', model_ws=model_ws,
verbose=True, exe_name='mt3dms')
mt2.change_model_ws(model_ws2)
mf2.write_input()
mt2.write_input()
success, buff = mf2.run_model(report=True)
assert success, 'MODFLOW did not run'
success, buff = mt2.run_model(report=True, normal_msg='Program completed.')
assert success, 'MT3D did not run'

fname = os.path.join(model_ws, 'MT3D001.UCN')
ucnobj = flopy.utils.UcnFile(fname)
conca = ucnobj.get_alldata()

fname = os.path.join(model_ws2, 'MT3D001.UCN')
ucnobj = flopy.utils.UcnFile(fname)
concb = ucnobj.get_alldata()

assert np.allclose(conca, concb)

return


def test_none_spdtype():
# ensure that -1 and None work as valid list entries in the
# stress period dictionary
model_ws = os.path.join('.', 'temp', 't068c')
mf = flopy.modflow.Modflow(model_ws=model_ws, exe_name=mf_exe_name)
dis = flopy.modflow.ModflowDis(mf, nper=2)
bas = flopy.modflow.ModflowBas(mf)
lpf = flopy.modflow.ModflowLpf(mf)
spd = {0: -1, 1: None}
wel = flopy.modflow.ModflowWel(mf, stress_period_data=spd)
pcg = flopy.modflow.ModflowPcg(mf)
mf.write_input()
mf2 = flopy.modflow.Modflow.load('modflowtest.nam', model_ws=model_ws,
verbose=True)
if run:
success, buff = mf.run_model(report=True)
assert success


if __name__ == '__main__':
test_mt3d_ssm_with_nodata_in_1st_sp()
test_none_spdtype()


12 changes: 8 additions & 4 deletions flopy/mt3d/mtssm.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -682,12 +682,13 @@ def load(f, model, nlay=None, nrow=None, ncol=None, nper=None,
print(' loading NSS...')
line = f.readline()
nss = int(line[0:10])
if model.verbose:
print(' NSS {}'.format(nss))

# Item D8: KSS, ISS, JSS, CSS, ITYPE, (CSSMS(n),n=1,NCOMP)
if model.verbose:
print(
' loading KSS, ISS, JSS, CSS, ITYPE, (CSSMS(n),n=1,NCOMP)...')
current = 0
print(' loading KSS, ISS, JSS, CSS, ITYPE, '
'(CSSMS(n),n=1,NCOMP)...')
if nss > 0:
current = np.empty((nss), dtype=dtype)
for ibnd in range(nss):
Expand All @@ -708,7 +709,10 @@ def load(f, model, nlay=None, nrow=None, ncol=None, nper=None,
current['i'] -= 1
current['j'] -= 1
current = current.view(np.recarray)
stress_period_data[iper] = current
stress_period_data[iper] = current
elif nss == 0:
stress_period_data[iper] = nss


if openfile:
f.close()
Expand Down
6 changes: 5 additions & 1 deletion flopy/pakbase.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -809,6 +809,7 @@ def load(f, model, pak_type, ext_unit_dict=None, **kwargs):
# read data for every stress period
bnd_output = None
stress_period_data = {}
current = None
for iper in range(nper):
if model.verbose:
msg = ' loading ' + str(pak_type) + \
Expand Down Expand Up @@ -843,7 +844,10 @@ def load(f, model, pak_type, ext_unit_dict=None, **kwargs):
current['node'] -= 1
bnd_output = np.recarray.copy(current)
else:
bnd_output = np.recarray.copy(current)
if current is None:
bnd_output = None
else:
bnd_output = np.recarray.copy(current)

for iparm in range(itmpp):
line = f.readline()
Expand Down
8 changes: 5 additions & 3 deletions flopy/utils/util_list.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -228,6 +228,8 @@ def dtype(self):
def get_itmp(self, kper):
if kper not in list(self.__data.keys()):
return None
if self.__vtype[kper] is None:
return -1
# If an external file, have to load it
if self.__vtype[kper] == str:
return self.__fromfile(self.__data[kper]).shape[0]
Expand Down Expand Up @@ -338,6 +340,9 @@ def __cast_data(self, data):
self.__cast_int(kper, d)
elif isinstance(d, str):
self.__cast_str(kper, d)
elif d is None:
self.__data[kper] = -1
self.__vtype[kper] = None
else:
raise Exception("MfList error: unsupported data type: " +
str(type(d)) + " at kper " +
Expand Down Expand Up @@ -374,9 +379,6 @@ def __cast_int(self, kper, d):
self.__data[kper] = 0
self.__vtype[kper] = None
else:
if kper == 0:
raise Exception("MfList error: dict integer value for " + \
"kper 0 for cannot be negative")
self.__data[kper] = -1
self.__vtype[kper] = None

Expand Down

0 comments on commit 51d8064

Please sign in to comment.