Skip to content

Commit

Permalink
update(remap_array): trap for None type idomain (#2034)
Browse files Browse the repository at this point in the history
* reset mfarray value for idomain to constant 1 if idomain is unspecified
  • Loading branch information
jlarsen-usgs authored Dec 7, 2023
1 parent 44abb51 commit 5a4533b
Show file tree
Hide file tree
Showing 2 changed files with 113 additions and 1 deletion.
109 changes: 109 additions & 0 deletions autotest/test_model_splitter.py
Original file line number Diff line number Diff line change
Expand Up @@ -582,3 +582,112 @@ def test_transient_array(function_tmpdir):
"storage package transient dictionary "
+ f"does not match for model '{name}'"
)


@requires_exe("mf6")
def test_idomain_none(function_tmpdir):
name = "id_test"
sim_path = function_tmpdir
new_sim_path = function_tmpdir / f"{name}_split_model"

tdis_data = [
(0., 1, 1.0),
(300000., 1, 1.0),
(36500., 10, 1.5),
(300000, 1, 1.0),
]

nper = len(tdis_data)
nlay, nrow, ncol = 3, 21, 20
xlen, ylen = 10000., 10500.
delc = ylen / nrow
delr = xlen / ncol

top = 400.
botm = [220., 200., 0]
K11 = [50., 0.01, 200.]
K33 = [10., 0.01, 20.]
Ss, Sy = 0.0001, 0.1
H_east = 320.
recharge = 0.005
idomain = None

sim = flopy.mf6.MFSimulation(sim_name=name, sim_ws=sim_path)
tdis = flopy.mf6.ModflowTdis(
sim, nper=nper, perioddata=tdis_data, time_units="days"
)
ims = flopy.mf6.ModflowIms(
sim,
complexity="simple",
print_option="all",
outer_dvclose=1e-6,
inner_dvclose=1e-6
)

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=top,
botm=botm,
length_units="meters",
idomain=idomain
)
npf = flopy.mf6.ModflowGwfnpf(
gwf, icelltype=1, save_specific_discharge=True, k=K11, k33=K33
)
sto = flopy.mf6.ModflowGwfsto(
gwf,
iconvert=1,
ss=Ss,
sy=Sy,
steady_state={0: True, 3: True},
transient={2: True}
)
ic = flopy.mf6.ModflowGwfic(gwf, strt=top)
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=recharge)

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

well_spd = {
1: [(0, 10, 9, -75000.)],
2: [
(0, 10, 9, -75000.),
(2, 12, 4, -100000.)
],
}
wel = flopy.mf6.ModflowGwfwel(gwf, stress_period_data=well_spd)
sim.write_simulation()
sim.run_simulation()

ms = Mf6Splitter(sim)
sarr = ms.optimize_splitting_mask(3)
new_sim = ms.split_model(sarr)
new_sim.set_sim_path(new_sim_path)
new_sim.write_simulation()
new_sim.run_simulation()

kstpkper = (9, 2)
head = gwf.output.head().get_data(kstpkper=kstpkper)
head_dict = {}
for idx, modelname in enumerate(new_sim.model_names):
mnum = int(modelname.split("_")[-1])
h = new_sim.get_model(modelname).output.head().get_data(
kstpkper=kstpkper)
head_dict[mnum] = h
new_head = ms.reconstruct_array(head_dict)

err_msg = "Heads from original and split models do not match"
np.testing.assert_allclose(new_head, head, atol=1e-07, err_msg=err_msg)

5 changes: 4 additions & 1 deletion flopy/mf6/utils/model_splitter.py
Original file line number Diff line number Diff line change
Expand Up @@ -1128,7 +1128,10 @@ def _remap_array(self, item, mfarray, mapped_data, **kwargs):
fnames = kwargs.pop("fnames", None)
if not hasattr(mfarray, "size"):
if mfarray.array is None:
return mapped_data
if item == "idomain":
mfarray.set_data(1)
else:
return mapped_data

how = [
i.data_storage_type.value
Expand Down

0 comments on commit 5a4533b

Please sign in to comment.