From 51d80641e81f844490306d2cddefacb21bda8302 Mon Sep 17 00:00:00 2001 From: langevin-usgs Date: Mon, 16 Dec 2019 15:17:19 -0600 Subject: [PATCH] fix(Mflist): allow None as a list entry (#765) * 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 --- autotest/t060_test_lkt.py | 71 ++++++++++++++++---------- autotest/t068_test_ssm.py | 103 +++++++++++++++++++++++++++++--------- flopy/mt3d/mtssm.py | 12 +++-- flopy/pakbase.py | 6 ++- flopy/utils/util_list.py | 8 +-- 5 files changed, 142 insertions(+), 58 deletions(-) diff --git a/autotest/t060_test_lkt.py b/autotest/t060_test_lkt.py index 0b077165f4..bd628bc4d6 100644 --- a/autotest/t060_test_lkt.py +++ b/autotest/t060_test_lkt.py @@ -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) @@ -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 @@ -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 @@ -124,16 +127,17 @@ 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 @@ -141,7 +145,7 @@ def test_lkt_with_multispecies(): # 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 @@ -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 @@ -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 @@ -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) @@ -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 @@ -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() diff --git a/autotest/t068_test_ssm.py b/autotest/t068_test_ssm.py index 45a6ca3060..89e34e8f5d 100644 --- a/autotest/t068_test_ssm.py +++ b/autotest/t068_test_ssm.py @@ -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 = {} @@ -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)) @@ -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) @@ -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() diff --git a/flopy/mt3d/mtssm.py b/flopy/mt3d/mtssm.py index a9451d46df..c7c7e7b3bb 100644 --- a/flopy/mt3d/mtssm.py +++ b/flopy/mt3d/mtssm.py @@ -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): @@ -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() diff --git a/flopy/pakbase.py b/flopy/pakbase.py index 15dfc052ec..95ef62b37d 100644 --- a/flopy/pakbase.py +++ b/flopy/pakbase.py @@ -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) + \ @@ -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() diff --git a/flopy/utils/util_list.py b/flopy/utils/util_list.py index 1abc0a4292..202d5d3d65 100644 --- a/flopy/utils/util_list.py +++ b/flopy/utils/util_list.py @@ -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] @@ -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 " + @@ -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