Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update(model_splitter): add timeseries support #2403

Merged
merged 5 commits into from
Jan 3, 2025
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
198 changes: 184 additions & 14 deletions autotest/test_model_splitter.py
Original file line number Diff line number Diff line change
Expand Up @@ -306,7 +306,9 @@ def test_control_records(function_tmpdir):
],
)

wel_rec = [((0, 4, 5), -10)]
wel_rec = [
((0, 4, 5), -10),
]

spd = {
0: wel_rec,
Expand Down Expand Up @@ -380,8 +382,22 @@ def test_empty_packages(function_tmpdir):
k33=20.0,
)
ic = flopy.mf6.ModflowGwfic(gwf, strt=0.0)
chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data={0: [((0, 0, 13), 0.0)]})
wel = flopy.mf6.ModflowGwfwel(gwf, stress_period_data={0: [((0, 0, 0), 1.0)]})
chd = flopy.mf6.ModflowGwfchd(
gwf,
stress_period_data={
0: [
((0, 0, 13), 0.0),
]
},
)
wel = flopy.mf6.ModflowGwfwel(
gwf,
stress_period_data={
0: [
((0, 0, 0), 1.0),
]
},
)

# Build SFR records
packagedata = [
Expand Down Expand Up @@ -429,7 +445,11 @@ def test_empty_packages(function_tmpdir):
nreaches=14,
packagedata=packagedata,
connectiondata=connectiondata,
perioddata={0: [(0, "INFLOW", 1.0)]},
perioddata={
0: [
(0, "INFLOW", 1.0),
]
},
)

array = np.zeros((nrow, ncol), dtype=int)
Expand Down Expand Up @@ -539,7 +559,10 @@ def test_transient_array(function_tmpdir):
for name in new_sim.model_names:
g = new_sim.get_model(name)
d = {}
for key in (0, 2):
for key in (
0,
2,
):
d[key] = g.sto.steady_state.get_data(key)
assert d == steady, (
"storage steady_state dictionary " + f"does not match for model '{name}'"
Expand Down Expand Up @@ -692,11 +715,39 @@ def test_unstructured_complex_disu(function_tmpdir):
iac, ja, ihc, cl12, hwva, angldegx = [], [], [], [], [], []
for cell, neigh in neighbors.items():
iac.append(len(neigh) + 1)
ihc.extend([1] * (len(neigh) + 1))
ja.extend([cell] + neigh)
cl12.extend([0] + [1] * len(neigh))
hwva.extend([0] + [1] * len(neigh))
adx = [0]
ihc.extend(
[
1,
]
* (len(neigh) + 1)
)
ja.extend(
[
cell,
]
+ neigh
)
cl12.extend(
[
0,
]
+ [
1,
]
* len(neigh)
)
hwva.extend(
[
0,
]
+ [
1,
]
* len(neigh)
)
adx = [
0,
]
for n in neigh:
ev = cell - n
if ev == -1 * ncol:
Expand Down Expand Up @@ -856,7 +907,12 @@ def string2geom(geostring, conversion=None):
)

ixs = flopy.utils.GridIntersect(modelgrid, method="vertex", rtree=True)
result = ixs.intersect([boundary], shapetype="Polygon")
result = ixs.intersect(
[
boundary,
],
shapetype="Polygon",
)
r, c = list(zip(*list(result.cellids)))
idomain = np.zeros(modelgrid.shape, dtype=int)
idomain[:, r, c] = 1
Expand Down Expand Up @@ -920,7 +976,12 @@ def string2geom(geostring, conversion=None):
botm[idx] = topc[idx] - dv0

strt = np.tile([modelgrid.top], (nlay, 1, 1))
idomain = np.tile([modelgrid.idomain[0]], (5, 1, 1))
idomain = np.tile(
[
modelgrid.idomain[0],
],
(5, 1, 1),
)

# setup recharge
dist_from_riv = 10000.0
Expand Down Expand Up @@ -1203,7 +1264,7 @@ def build_gwt_model(sim, gwtname, rch_package):
new_sim.run_simulation()

# compare results for each of the models
splits = list(range(nparts))
splits = range(nparts)
for name in sim.model_names:
gwm = sim.get_model(name)
if "concentration()" in gwm.output.methods():
Expand Down Expand Up @@ -1234,6 +1295,115 @@ def build_gwt_model(sim, gwtname, rch_package):
diff = np.nansum(diff)
if diff > 10.25:
raise AssertionError(
"Difference between output arrays: "
f"Difference between output arrays: "
f"{diff :.2f} greater than tolerance"
)


@requires_exe("mf6")
def test_timeseries(function_tmpdir):
sim = MFSimulation(
sim_name="np001",
sim_ws=function_tmpdir,
continue_=True,
memory_print_option="summary",
)

tdis_rc = [(6.0, 2, 1.0), (6.0, 3, 1.0)]
tdis = flopy.mf6.ModflowTdis(sim, time_units="DAYS", nper=2, perioddata=tdis_rc)

ims = flopy.mf6.ModflowIms(
sim,
print_option="ALL",
complexity="SIMPLE",
outer_dvclose=0.00001,
outer_maximum=50,
under_relaxation="NONE",
inner_maximum=30,
inner_dvclose=0.00001,
linear_acceleration="CG",
preconditioner_levels=7,
preconditioner_drop_tolerance=0.01,
number_orthogonalizations=2,
)
gwf = flopy.mf6.ModflowGwf(
sim,
)
dis = flopy.mf6.ModflowGwfdis(
gwf, nlay=1, nrow=1, ncol=10, delr=500, delc=500, top=100, botm=50
)

ic = flopy.mf6.ModflowGwfic(gwf, strt=95)

npf = flopy.mf6.ModflowGwfnpf(
gwf,
pname="npf_1",
save_flows=True,
alternative_cell_averaging="logarithmic",
icelltype=1,
k=5.0,
)

oc = flopy.mf6.ModflowGwfoc(
gwf,
budget_filerecord=[("np001_mod 1.cbc",)],
head_filerecord=[("np001_mod 1.hds",)],
saverecord={
0: [("HEAD", "ALL"), ("BUDGET", "ALL")],
1: [("HEAD", "ALL"), ("BUDGET", "ALL")],
},
printrecord=[("HEAD", "ALL")],
)

sto = flopy.mf6.ModflowGwfsto(
gwf, save_flows=True, iconvert=1, ss=0.000001, sy=0.15
)

wel = flopy.mf6.ModflowGwfwel(
gwf,
print_input=True,
print_flows=True,
save_flows=True,
maxbound=2,
stress_period_data={0: [(0, 0, 4, -2000.0), (0, 0, 7, -2.0)], 1: None},
)

tsdict = {
"filename": "drn_ts.ts",
"timeseries": [(0.0, 60.0), (100000.0, 60.0)],
"time_series_namerecord": "drn_1",
"interpolation_methodrecord": "linearend",
}
drn = flopy.mf6.ModflowGwfdrn(
gwf,
print_input=True,
print_flows=True,
save_flows=True,
maxbound=1,
timeseries=tsdict,
stress_period_data=[((0, 0, 0), 80, "drn_1")],
)

spd = {0: [((0, 0, 9), 110, 90.0, 100.0, 1.0, 2.0, 3.0)]}
riv = flopy.mf6.ModflowGwfriv(
gwf,
print_input=True,
print_flows=True,
save_flows=True,
maxbound=1,
auxiliary=["var1", "var2", "var3"],
stress_period_data=spd,
)
sim.write_simulation()
sim.run_simulation()

mfs = Mf6Splitter(sim)
mask = mfs.optimize_splitting_mask(2)
new_sim = mfs.split_model(mask)

new_sim.set_sim_path(function_tmpdir / "split_model")
new_sim.write_simulation()
success, _ = new_sim.run_simulation()

if not success:
raise AssertionError("Timeseries split simulation did not properly run")
68 changes: 57 additions & 11 deletions flopy/mf6/utils/model_splitter.py
Original file line number Diff line number Diff line change
Expand Up @@ -1025,6 +1025,7 @@ def _remap_filerecords(self, item, value, mapped_data):
"stage_filerecord",
"obs_filerecord",
"concentration_filerecord",
"ts_filerecord"
):
value = value.array
if value is None:
Expand Down Expand Up @@ -2111,6 +2112,33 @@ def _remap_csub(self, package, mapped_data):

return mapped_data

def _remap_buy(self, package, mapped_data):
"""
Method to remap a Buoyancy package

Parameters
----------
package : ModflowGwfbuy
mapped_data : dict
dictionary of remapped package data

Returns
-------
dict
"""
recarray = package.packagedata.array
mnames = [i for i in recarray["modelname"]]

for mkey in mapped_data.keys():
new_recarray = recarray.copy()
new_mnames = [
f"{mn}_{mkey :0{self._fdigits}d}" for mn in mnames
]
new_recarray["modelname"] = new_mnames
mapped_data[mkey]["packagedata"] = new_recarray

return mapped_data

def _set_boundname_remaps(self, recarray, obs_map, variables, mkey):
"""

Expand Down Expand Up @@ -2235,18 +2263,19 @@ def _remap_fmi(self, package, mapped_data):
dict
"""
packagedata = package.packagedata.array
fnames = packagedata.fname
if packagedata is not None:
fnames = packagedata.fname

for mkey in self._model_dict.keys():
new_fnames = []
for fname in fnames:
new_val = fname.split(".")
new_val = f"{'.'.join(new_val[0:-1])}_{mkey :0{self._fdigits}d}.{new_val[-1]}"
new_fnames.append(new_val)

new_packagedata = packagedata.copy()
new_packagedata["fname"] = new_fnames
mapped_data[mkey]["packagedata"] = new_packagedata
for mkey in self._model_dict.keys():
new_fnames = []
for fname in fnames:
new_val = fname.split(".")
new_val = f"{'.'.join(new_val[0:-1])}_{mkey :0{self._fdigits}d}.{new_val[-1]}"
new_fnames.append(new_val)

new_packagedata = packagedata.copy()
new_packagedata["fname"] = new_fnames
mapped_data[mkey]["packagedata"] = new_packagedata

return mapped_data

Expand Down Expand Up @@ -2888,6 +2917,9 @@ def _remap_package(self, package, ismvr=False):
elif isinstance(package, modflow.ModflowGwfcsub):
mapped_data = self._remap_csub(package, mapped_data)

elif isinstance(package, modflow.ModflowGwfbuy):
mapped_data = self._remap_buy(package, mapped_data)

elif isinstance(
package, (modflow.ModflowGwfuzf, modflow.ModflowGwtuzt)
):
Expand Down Expand Up @@ -3025,6 +3057,20 @@ def _remap_package(self, package, ismvr=False):
elif isinstance(value, mfdatascalar.MFScalar):
for mkey in self._model_dict.keys():
mapped_data[mkey][item] = value.data
elif isinstance(value, modflow.mfutlts.UtltsPackages):
if value._filerecord.array is None:
continue
tspkg = value._packages[0]
for mkey in self._model_dict.keys():
new_fname = tspkg.filename.split(".")
new_fname = f"{'.'.join(new_fname[0:-1])}_{mkey :0{self._fdigits}d}.{new_fname[-1]}"
tsdict = {
"filename": new_fname,
"timeseries": tspkg.timeseries.array,
"time_series_namerecord": tspkg.time_series_namerecord.array["time_series_names"][0],
"interpolation_methodrecord": tspkg.interpolation_methodrecord.array["interpolation_method"][0]
}
mapped_data[mkey]["timeseries"] = tsdict
else:
pass

Expand Down
Loading