Skip to content

Commit

Permalink
feat(apt): generalization of advanced package transport (#327)
Browse files Browse the repository at this point in the history 
* fix(uzf): fix indexing error in UZF (#274)

* Error introduced as part of recent UZF refactoring
* Closes #273

* fix(memory): some variables not deallocated (#278)

Implemented new check in develop mode so code bombs with error if memory manager variable not deallocated

* refactor(budobj): new budget object for advanced packages (#279)

* single code base for writing binary budget files for advanced packages
* single code base for creating and writing budget tables to list file for advanced packages
* implemented for MAW, UZF, LAK, SFR, and MVR
* closes #277
* update mf6exes from 2.0 to 3.0
* will allow generalized transport calculations for advanced packages

* refactor(advanced packages): read static data as part of df() (#281)

* refactor(advanced packages): modify advanced packages to read all static data as part of df()
* modify setup_budobj to include the connectivity so that it is available to other models

* fix(budterm): initialize nlist to zero (#283)

* initialize nlist to zero (#284)

* feat(mwt): generalization of lkt with more to do

* update mwt comment table

* working to generalize the gwt1apt1, should still fail for mwt problem

* heavy changes to get an advanced package transport module working

* more generalization for apt; working for mvt now

* cleanup on apt
  • Loading branch information
@langevin-usgs
langevin-usgs committed Mar 6, 2020
1 parent f55cf79 commit 5d8d57c
Show file tree
Hide file tree
Showing 25 changed files with 5,615 additions and 2,885 deletions.
390 changes: 390 additions & 0 deletions autotest/test_gwt_lkt03.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,390 @@
# Simple one-layer model with a lak. Purpose is to test the lak flow terms
# such as rainfall to make sure mixing is calculated correctly.

import os
import sys
import numpy as np

try:
import flopy
except:
msg = 'Error. FloPy package is not available.\n'
msg += 'Try installing using the following command:\n'
msg += ' pip install flopy'
raise Exception(msg)

from framework import testing_framework
from simulation import Simulation

ex = ['lkt_03']
exdirs = []
for s in ex:
exdirs.append(os.path.join('temp', s))


def get_model(idx, dir):
lx = 7.
lz = 1.
nlay = 1
nrow = 1
ncol = 7
nper = 1
delc = 1.
delr = lx / ncol
delz = lz / nlay
top = [0., 0., -0.10, -0.10, -0.10, 0., 0.]
botm = [[-1, -1, -1, -1, -1, -1, -1]]

perlen = [1.0]
nstp = [10]
tsmult = [1.]

Kh = 20.
Kv = 20.

steady = [True]
tdis_rc = []
for i in range(nper):
tdis_rc.append((perlen[i], nstp[i], tsmult[i]))

single_matrix = False
nouter, ninner = 700, 300
hclose, rclose, relax = 1e-8, 1e-6, 0.97

name = ex[idx]

# build MODFLOW 6 files
ws = dir
sim = flopy.mf6.MFSimulation(sim_name=name, version='mf6',
exe_name='mf6',
sim_ws=ws)
# create tdis package
tdis = flopy.mf6.ModflowTdis(sim, time_units='DAYS',
nper=nper, perioddata=tdis_rc)

# create gwf model
gwfname = 'gwf_' + name

gwf = flopy.mf6.MFModel(sim, model_type='gwf6', modelname=gwfname,
model_nam_file='{}.nam'.format(gwfname))

imsgwf = flopy.mf6.ModflowIms(sim, print_option='ALL',
outer_hclose=hclose,
outer_maximum=nouter,
under_relaxation='NONE',
inner_maximum=ninner,
inner_hclose=hclose, rcloserecord=rclose,
linear_acceleration='BICGSTAB',
scaling_method='NONE',
reordering_method='NONE',
relaxation_factor=relax,
filename='{}.ims'.format(gwfname))

idomain = np.full((nlay, nrow, ncol), 1)
dis = flopy.mf6.ModflowGwfdis(gwf, nlay=nlay, nrow=nrow, ncol=ncol,
delr=delr, delc=delc,
top=top, botm=botm, idomain=idomain)

# initial conditions
ic = flopy.mf6.ModflowGwfic(gwf, strt=0.)

# node property flow
npf = flopy.mf6.ModflowGwfnpf(gwf, xt3doptions=False,
save_flows=True,
save_specific_discharge=True,
icelltype=0,
k=Kh, k33=Kv)

# chd files
chdlist1 = [[(0, 0, 0), 0.0, 0.],
[(0, 0, ncol - 1), 0., 0.],
]
chd1 = flopy.mf6.ModflowGwfchd(gwf,
stress_period_data=chdlist1,
print_input=True,
print_flows=True,
save_flows=False,
pname='CHD-1',
auxiliary='CONCENTRATION',
filename='{}.chd'.format(gwfname))

# pak_data = [lakeno, strt, nlakeconn, CONC, dense, boundname]
pak_data = [(0, 0., 1, 0., 1025.),
(1, 0., 1, 0., 1025.),
(2, 0., 1, 0., 1025.)]

connlen = connwidth = delr / 2.
con_data = []
# con_data=(lakeno,iconn,(cellid),claktype,bedleak,belev,telev,connlen,connwidth )
# lake 1
con_data.append(
(0, 0, (0, 0, 2), 'VERTICAL', 'None', 10, 10, connlen, connwidth))
# lake 2
con_data.append(
(1, 0, (0, 0, 3), 'VERTICAL', 'None', 10, 10, connlen, connwidth))
# lake 3
con_data.append(
(2, 0, (0, 0, 4), 'VERTICAL', 'None', 10, 10, connlen, connwidth))

p_data = [
(0, 'RAINFALL', .1),
(1, 'RAINFALL', .1),
(2, 'RAINFALL', .1),
(0, 'STATUS', 'CONSTANT'),
(1, 'STATUS', 'CONSTANT'),
(2, 'STATUS', 'CONSTANT'),
]

# note: for specifying lake number, use fortran indexing!
lak_obs = {('lak_obs.csv'): [('lakestage', 'stage', 1),
('lakevolume', 'volume', 1),
('lak1', 'lak', 1, 1),
('lak2', 'lak', 1, 2),
('lak3', 'lak', 1, 3)]}

lak = flopy.mf6.modflow.ModflowGwflak(gwf, save_flows=True,
print_input=True,
print_flows=True,
print_stage=True,
stage_filerecord=gwfname+'.lak.stg',
budget_filerecord=gwfname+'.lak.bud',
nlakes=3, ntables=0, noutlets=0,
packagedata=pak_data,
#outlets=outlets,
pname='LAK-1',
connectiondata=con_data,
lakeperioddata=p_data,
#outletperioddata=outletperioddata,
observations=lak_obs,
auxiliary=['CONCENTRATION',
'DENSITY'])

# output control
oc = flopy.mf6.ModflowGwfoc(gwf,
budget_filerecord='{}.cbc'.format(gwfname),
head_filerecord='{}.hds'.format(gwfname),
headprintrecord=[
('COLUMNS', 10, 'WIDTH', 15,
'DIGITS', 6, 'GENERAL')],
saverecord=[('HEAD', 'ALL', 'STEPS'),
('BUDGET', 'ALL', 'STEPS')],
printrecord=[('HEAD', 'LAST', 'STEPS'),
('BUDGET', 'LAST', 'STEPS')])

# create gwt model
if True:
gwtname = 'gwt_' + name
gwt = flopy.mf6.MFModel(sim, model_type='gwt6', modelname=gwtname,
model_nam_file='{}.nam'.format(gwtname))

if not single_matrix:
imsgwt = flopy.mf6.ModflowIms(sim, print_option='ALL',
outer_hclose=hclose,
outer_maximum=nouter,
under_relaxation='NONE',
inner_maximum=ninner,
inner_hclose=hclose, rcloserecord=rclose,
linear_acceleration='BICGSTAB',
scaling_method='NONE',
reordering_method='NONE',
relaxation_factor=relax,
filename='{}.ims'.format(gwtname))
sim.register_ims_package(imsgwt, [gwt.name])

dis = flopy.mf6.ModflowGwtdis(gwt, nlay=nlay, nrow=nrow, ncol=ncol,
delr=delr, delc=delc,
top=top, botm=botm, idomain=idomain)

# initial conditions
ic = flopy.mf6.ModflowGwtic(gwt, strt=0.)

# advection
adv = flopy.mf6.ModflowGwtadv(gwt, scheme='UPSTREAM')

# dispersion
# diffc = 0.0
# dsp = flopy.mf6.ModflowGwtdsp(gwt, xt3d=True, diffc=diffc,
# alh=0.1, ath1=0.01, atv=0.05,
# filename='{}.dsp'.format(gwtname))

# storage
porosity = 0.30
sto = flopy.mf6.ModflowGwtmst(gwt, porosity=porosity)
# sources
sourcerecarray = [('CHD-1', 'AUX', 'CONCENTRATION'),
#('WEL-1', 'AUX', 'CONCENTRATION'),
]
ssm = flopy.mf6.ModflowGwtssm(gwt, sources=sourcerecarray,
filename='{}.ssm'.format(gwtname))

lktpackagedata = [(0, 0., 99., 999., 'mylake1'),
(1, 0., 99., 999., 'mylake2'),
(2, 0., 99., 999., 'mylake3'),]
lktperioddata = [
(0, 'STATUS', 'ACTIVE'),
(1, 'STATUS', 'ACTIVE'),
(2, 'STATUS', 'ACTIVE'),
(0, 'RAINFALL', 0.),
(1, 'RAINFALL', 1.),
(2, 'RAINFALL', 10.),
]

lkt_obs = {(gwtname + '.lkt.obs.csv', ):
[
('lkt-1-conc', 'CONCENTRATION', 1),
('lkt-1-extinflow', 'EXT-INFLOW', 1),
('lkt-1-rain', 'RAINFALL', 1),
('lkt-1-roff', 'RUNOFF', 1),
('lkt-1-evap', 'EVAPORATION', 1),
('lkt-1-wdrl', 'WITHDRAWAL', 1),
('lkt-1-stor', 'STORAGE', 1),
('lkt-1-const', 'CONSTANT', 1),
('lkt-1-gwt1', 'LKT', 1, 1),
('lkt-1-gwt2', 'LKT', 1, 2),
('lkt-2-gwt1', 'LKT', 2, 1),
('lkt-1-mylake1', 'LKT', 'MYLAKE1'),
('lkt-1-fjf', 'FLOW-JA-FACE', 1, 2),
('lkt-2-fjf', 'FLOW-JA-FACE', 2, 1),
('lkt-3-fjf', 'FLOW-JA-FACE', 2, 3),
('lkt-4-fjf', 'FLOW-JA-FACE', 3, 2),
('lkt-5-fjf', 'FLOW-JA-FACE', 'MYLAKE1'),
('lkt-6-fjf', 'FLOW-JA-FACE', 'MYLAKE2'),
('lkt-7-fjf', 'FLOW-JA-FACE', 'MYLAKE3'),
],
}
# append additional obs attributes to obs dictionary
lkt_obs['digits'] = 7
lkt_obs['print_input'] = True
lkt_obs['filename'] = gwtname + '.lkt.obs'

lkt = flopy.mf6.modflow.ModflowGwtlkt(gwt,
boundnames=True,
save_flows=True,
print_input=True,
print_flows=True,
print_concentration=True,
concentration_filerecord=gwtname + '.lkt.bin',
budget_filerecord=gwtname + '.lkt.bud',
packagedata=lktpackagedata,
lakeperioddata=lktperioddata,
#observations=lkt_obs,
pname='LAK-1',
auxiliary=['aux1', 'aux2'])
# output control
oc = flopy.mf6.ModflowGwtoc(gwt,
budget_filerecord='{}.cbc'.format(gwtname),
concentration_filerecord='{}.ucn'.format(
gwtname),
concentrationprintrecord=[
('COLUMNS', 10, 'WIDTH', 15,
'DIGITS', 6, 'GENERAL')],
saverecord=[('CONCENTRATION', 'ALL'),
('BUDGET', 'ALL')],
printrecord=[('CONCENTRATION', 'ALL'),
('BUDGET', 'ALL')])

# GWF GWT exchange
gwfgwt = flopy.mf6.ModflowGwfgwt(sim, exgtype='GWF6-GWT6',
exgmnamea=gwfname, exgmnameb=gwtname,
filename='{}.gwfgwt'.format(name))



return sim


def build_models():
for idx, dir in enumerate(exdirs):
sim = get_model(idx, dir)
sim.write_simulation()
return


def eval_results(sim):
print('evaluating results...')

# ensure lake concentrations were saved
name = ex[sim.idxsim]
gwtname = 'gwt_' + name
fname = gwtname + '.lkt.bin'
fname = os.path.join(sim.simpath, fname)
assert os.path.isfile(fname)

# load the lake concentrations and make sure all values are correct
cobj = flopy.utils.HeadFile(fname, text='CONCENTRATION')
clak = cobj.get_data()
answer = np.array([0., 1., 10.])
assert np.allclose(clak, answer), '{} {}'.format(clak, answer)

# load the aquifer concentrations and make sure all values are correct
fname = gwtname + '.ucn'
fname = os.path.join(sim.simpath, fname)
cobj = flopy.utils.HeadFile(fname, text='CONCENTRATION')
caq = cobj.get_data()
answer = np.zeros((7,))
assert np.allclose(caq, answer), '{} {}'.format(caq.flatten(), answer)

# load the lake budget file
fname = gwtname + '.lkt.bud'
fname = os.path.join(sim.simpath, fname)
assert os.path.isfile(fname)
bobj = flopy.utils.CellBudgetFile(fname, precision='double', verbose=False)

# check the rainfall terms
res = bobj.get_data(text='rainfall')[-1]
answer = [(1, 1, 0.), (2, 2, 0.1), (3, 3, 1.)]
dt = [('node', '<i4'), ('node2', '<i4'), ('q', '<f8')]
answer = np.array(answer, dtype=dt)
for dtname, dttype in dt:
assert np.allclose(res[dtname], answer[dtname]), '{} {}'.format(res, answer)

# check the storage terms, which include the total mass in the lake as an aux variable
res = bobj.get_data(text='storage')[-1]
answer = [(1, 1, 0., 0.), (2, 2, -0.1, 0.1), (3, 3, -1., 1.)]
dt = [('node', '<i4'), ('node2', '<i4'), ('q', '<f8'), ('MASS', '<f8')]
answer = np.array(answer, dtype=dt)
for dtname, dttype in dt:
assert np.allclose(res[dtname], answer[dtname]), '{} {}'.format(res, answer)

# uncomment when testing
# assert False

return


# - No need to change any code below
def test_mf6model():
# initialize testing framework
test = testing_framework()

# build the models
build_models()

# run the test models
for idx, dir in enumerate(exdirs):
yield test.run_mf6, Simulation(dir, exfunc=eval_results, idxsim=idx)

return


def main():
# initialize testing framework
test = testing_framework()

# build the models
build_models()

# run the test models
for idx, dir in enumerate(exdirs):
sim = Simulation(dir, exfunc=eval_results, idxsim=idx)
test.run_mf6(sim)

return


if __name__ == "__main__":
# print message
print('standalone run of {}'.format(os.path.basename(__file__)))

# run main routine
main()
Loading

0 comments on commit 5d8d57c

Please sign in to comment.