Skip to content

Commit

Permalink
fix(model_splitter): check keys in mftransient array (#1998)
Browse files Browse the repository at this point in the history 
* fix issue with MFScalarTransient data in splitter (storage package period data)
* add test for MFScalarTransient issue
  • Loading branch information
@jdhughes-usgs
jdhughes-usgs authored Nov 13, 2023
1 parent 696a209 commit 5ebc216
Show file tree
Hide file tree
Showing 2 changed files with 164 additions and 41 deletions.
168 changes: 139 additions & 29 deletions autotest/test_model_splitter.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -356,7 +356,7 @@ def test_control_records(function_tmpdir):

spd_ls1 = ml1.wel.stress_period_data.get_record(1)
spd_ls2 = ml1.wel.stress_period_data.get_record(2)

if spd_ls1["filename"] is None or spd_ls1["binary"]:
raise AssertionError(
"External ascii files not being preserved for MFList"
Expand All @@ -379,41 +379,49 @@ def test_empty_packages(function_tmpdir):
nrow, ncol = 1, 14
base_name = "sfr01gwfgwf"
gwf = flopy.mf6.ModflowGwf(sim, modelname=base_name, save_flows=True)
dis = flopy.mf6.ModflowGwfdis(gwf, nrow=1, ncol=14, top=0., botm=-1.)
dis = flopy.mf6.ModflowGwfdis(gwf, nrow=1, ncol=14, top=0.0, botm=-1.0)
npf = flopy.mf6.ModflowGwfnpf(
gwf,
save_flows=True,
save_specific_discharge=True,
icelltype=0,
k=20.,
k33=20.
k=20.0,
k33=20.0,
)
ic = flopy.mf6.ModflowGwfic(gwf, strt=0.)
ic = flopy.mf6.ModflowGwfic(gwf, strt=0.0)
chd = flopy.mf6.ModflowGwfchd(
gwf,
stress_period_data={0: [((0, 0, 13), 0.0),]}
stress_period_data={
0: [
((0, 0, 13), 0.0),
]
},
)
wel = flopy.mf6.ModflowGwfwel(
gwf,
stress_period_data={0: [((0, 0, 0), 1.),]}
stress_period_data={
0: [
((0, 0, 0), 1.0),
]
},
)

# Build SFR records
packagedata = [
(0, (0, 0, 0), 1., 1., 1., 0., 0.1, 0., 1., 1, 1., 0),
(1, (0, 0, 1), 1., 1., 1., 0., 0.1, 0., 1., 2, 1., 0),
(2, (0, 0, 2), 1., 1., 1., 0., 0.1, 0., 1., 2, 1., 0),
(3, (0, 0, 3), 1., 1., 1., 0., 0.1, 0., 1., 2, 1., 0),
(4, (0, 0, 4), 1., 1., 1., 0., 0.1, 0., 1., 2, 1., 0),
(5, (0, 0, 5), 1., 1., 1., 0., 0.1, 0., 1., 2, 1., 0),
(6, (0, 0, 6), 1., 1., 1., 0., 0.1, 0., 1., 2, 1., 0),
(7, (0, 0, 7), 1., 1., 1., 0., 0.1, 0., 1., 2, 1., 0),
(8, (0, 0, 8), 1., 1., 1., 0., 0.1, 0., 1., 2, 1., 0),
(9, (0, 0, 9), 1., 1., 1., 0., 0.1, 0., 1., 2, 1., 0),
(10, (0, 0, 10), 1., 1., 1., 0., 0.1, 0., 1., 2, 1., 0),
(11, (0, 0, 11), 1., 1., 1., 0., 0.1, 0., 1., 2, 1., 0),
(12, (0, 0, 12), 1., 1., 1., 0., 0.1, 0., 1., 2, 1., 0),
(13, (0, 0, 13), 1., 1., 1., 0., 0.1, 0., 1., 1, 1., 0),
(0, (0, 0, 0), 1.0, 1.0, 1.0, 0.0, 0.1, 0.0, 1.0, 1, 1.0, 0),
(1, (0, 0, 1), 1.0, 1.0, 1.0, 0.0, 0.1, 0.0, 1.0, 2, 1.0, 0),
(2, (0, 0, 2), 1.0, 1.0, 1.0, 0.0, 0.1, 0.0, 1.0, 2, 1.0, 0),
(3, (0, 0, 3), 1.0, 1.0, 1.0, 0.0, 0.1, 0.0, 1.0, 2, 1.0, 0),
(4, (0, 0, 4), 1.0, 1.0, 1.0, 0.0, 0.1, 0.0, 1.0, 2, 1.0, 0),
(5, (0, 0, 5), 1.0, 1.0, 1.0, 0.0, 0.1, 0.0, 1.0, 2, 1.0, 0),
(6, (0, 0, 6), 1.0, 1.0, 1.0, 0.0, 0.1, 0.0, 1.0, 2, 1.0, 0),
(7, (0, 0, 7), 1.0, 1.0, 1.0, 0.0, 0.1, 0.0, 1.0, 2, 1.0, 0),
(8, (0, 0, 8), 1.0, 1.0, 1.0, 0.0, 0.1, 0.0, 1.0, 2, 1.0, 0),
(9, (0, 0, 9), 1.0, 1.0, 1.0, 0.0, 0.1, 0.0, 1.0, 2, 1.0, 0),
(10, (0, 0, 10), 1.0, 1.0, 1.0, 0.0, 0.1, 0.0, 1.0, 2, 1.0, 0),
(11, (0, 0, 11), 1.0, 1.0, 1.0, 0.0, 0.1, 0.0, 1.0, 2, 1.0, 0),
(12, (0, 0, 12), 1.0, 1.0, 1.0, 0.0, 0.1, 0.0, 1.0, 2, 1.0, 0),
(13, (0, 0, 13), 1.0, 1.0, 1.0, 0.0, 0.1, 0.0, 1.0, 1, 1.0, 0),
]

connectiondata = [
Expand All @@ -430,7 +438,7 @@ def test_empty_packages(function_tmpdir):
(10, 9, -11),
(11, 10, -12),
(12, 11, -13),
(13, 12)
(13, 12),
]

sfr = flopy.mf6.ModflowGwfsfr(
Expand All @@ -444,7 +452,11 @@ def test_empty_packages(function_tmpdir):
nreaches=14,
packagedata=packagedata,
connectiondata=connectiondata,
perioddata={0: [(0, "INFLOW", 1.),]}
perioddata={
0: [
(0, "INFLOW", 1.0),
]
},
)

array = np.zeros((nrow, ncol), dtype=int)
Expand All @@ -456,14 +468,10 @@ def test_empty_packages(function_tmpdir):
m1 = new_sim.get_model(f"{base_name}_1")

if "chd_0" in m0.package_dict:
raise AssertionError(
f"Empty CHD file written to {base_name}_0 model"
)
raise AssertionError(f"Empty CHD file written to {base_name}_0 model")

if "wel_0" in m1.package_dict:
raise AssertionError(
f"Empty WEL file written to {base_name}_1 model"
)
raise AssertionError(f"Empty WEL file written to {base_name}_1 model")

mvr_status0 = m0.sfr.mover.array
mvr_status1 = m0.sfr.mover.array
Expand All @@ -472,3 +480,105 @@ def test_empty_packages(function_tmpdir):
raise AssertionError(
"Mover status being overwritten in options splitting"
)


@requires_exe("mf6")
def test_transient_array(function_tmpdir):
name = "tarr"
new_sim_path = function_tmpdir / f"{name}_split_model"
nper = 3
tdis_data = [
(300000.0, 1, 1.0),
(36500.0, 10, 1.5),
(300000, 1, 1.0),
]
nlay, nrow, ncol = 3, 21, 20
xlen, ylen = 10000.0, 10500.0
delc = ylen / nrow
delr = xlen / ncol
steady = {0: True, 2: True}
transient = {1: True}

sim = flopy.mf6.MFSimulation(sim_name=name, sim_ws=new_sim_path)
ims = flopy.mf6.ModflowIms(sim, complexity="simple")
tdis = flopy.mf6.ModflowTdis(sim, nper=nper, perioddata=tdis_data)

gwf = flopy.mf6.ModflowGwf(
sim,
modelname=name,
save_flows=True,
)
dis = flopy.mf6.ModflowGwfdis(
gwf,
nlay=nlay,
nrow=nrow,
ncol=ncol,
delr=delr,
delc=delc,
top=400,
botm=[220.0, 200.0, 0],
length_units="meters",
)
npf = flopy.mf6.ModflowGwfnpf(
gwf,
icelltype=1,
k=[50.0, 0.01, 200.0],
k33=[10.0, 0.01, 20.0],
)
sto = flopy.mf6.ModflowGwfsto(
gwf,
iconvert=1,
ss=0.0001,
sy=0.1,
steady_state=steady,
transient=transient,
)
ic = flopy.mf6.ModflowGwfic(gwf, strt=400.0)
oc = flopy.mf6.ModflowGwfoc(
gwf,
head_filerecord=f"{name}.hds",
budget_filerecord=f"{name}.cbc",
saverecord={0: [("head", "all"), ("budget", "all")]},
)

rch = flopy.mf6.ModflowGwfrcha(gwf, recharge=0.005)

chd_spd = [(0, i, ncol - 1, 320) for i in range(nrow)]
chd = flopy.mf6.ModflowGwfchd(
gwf,
stress_period_data=chd_spd,
)

well_spd = {
0: [(0, 10, 9, -75000.0)],
2: [(0, 10, 9, -75000.0), (2, 12, 4, -100000.0)],
}
wel = flopy.mf6.ModflowGwfwel(
gwf,
stress_period_data=well_spd,
)

sarr = np.ones((nrow, ncol), dtype=int)
sarr[:, int(ncol / 2) :] = 2
mfsplit = Mf6Splitter(sim)
new_sim = mfsplit.split_model(sarr)

for name in new_sim.model_names:
g = new_sim.get_model(name)
d = {}
for key in (
0,
2,
):
d[key] = g.sto.steady_state.get_data(key).get_data()
assert d == steady, (
"storage steady_state dictionary "
+ f"does not match for model '{name}'"
)
d = {}
for key in (1,):
d[key] = g.sto.transient.get_data(key).get_data()
assert d == transient, (
"storage package transient dictionary "
+ f"does not match for model '{name}'"
)
37 changes: 25 additions & 12 deletions flopy/mf6/utils/model_splitter.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1076,20 +1076,30 @@ def _remap_transient_array(self, item, mftransient, mapped_data):

d0 = {mkey: {} for mkey in self._model_dict.keys()}
for per, array in enumerate(mftransient.array):
storage = mftransient._data_storage[per]
how = [
i.data_storage_type.value
for i in storage.layer_storage.multi_dim_list
]
binary = [i.binary for i in storage.layer_storage.multi_dim_list]
fnames = [i.fname for i in storage.layer_storage.multi_dim_list]
if per in mftransient._data_storage.keys():
storage = mftransient._data_storage[per]
how = [
i.data_storage_type.value
for i in storage.layer_storage.multi_dim_list
]
binary = [
i.binary for i in storage.layer_storage.multi_dim_list
]
fnames = [
i.fname for i in storage.layer_storage.multi_dim_list
]

d = self._remap_array(
item, array, mapped_data, how=how, binary=binary, fnames=fnames
)
d = self._remap_array(
item,
array,
mapped_data,
how=how,
binary=binary,
fnames=fnames,
)

for mkey in d.keys():
d0[mkey][per] = d[mkey][item]
for mkey in d.keys():
d0[mkey][per] = d[mkey][item]

for mkey, values in d0.items():
mapped_data[mkey][item] = values
Expand Down Expand Up @@ -2936,6 +2946,9 @@ def _remap_package(self, package, ismvr=False):
value = list_data
mapped_data = self._remap_mflist(item, value, mapped_data)

elif isinstance(value, mfdatascalar.MFScalarTransient):
for mkey in self._model_dict.keys():
mapped_data[mkey][item] = value._data_storage
elif isinstance(value, mfdatascalar.MFScalar):
for mkey in self._model_dict.keys():
mapped_data[mkey][item] = value.data
Expand Down

0 comments on commit 5ebc216

Please sign in to comment.