From e5bd078765e02e4a272ebf69d6a8917baf22fced Mon Sep 17 00:00:00 2001 From: jlarsen Date: Fri, 20 Dec 2024 15:39:05 -0800 Subject: [PATCH 1/5] Update(model_splitter): add timeseries support * add support for empty fmi package * add support for buy package mapping --- flopy/mf6/utils/model_splitter.py | 68 ++++++++++++++++++++++++++----- 1 file changed, 57 insertions(+), 11 deletions(-) diff --git a/flopy/mf6/utils/model_splitter.py b/flopy/mf6/utils/model_splitter.py index a4db88e75..9022d9844 100644 --- a/flopy/mf6/utils/model_splitter.py +++ b/flopy/mf6/utils/model_splitter.py @@ -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: @@ -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): """ @@ -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 @@ -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) ): @@ -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 From cee55119e4ddb5e07feb37fcc81f7ef0e14fa1a1 Mon Sep 17 00:00:00 2001 From: jlarsen Date: Fri, 20 Dec 2024 15:39:35 -0800 Subject: [PATCH 2/5] add autotest --- autotest/test_model_splitter.py | 206 +++++++++++++++++++++++++++++--- 1 file changed, 191 insertions(+), 15 deletions(-) diff --git a/autotest/test_model_splitter.py b/autotest/test_model_splitter.py index c2ddeab87..590c6f894 100644 --- a/autotest/test_model_splitter.py +++ b/autotest/test_model_splitter.py @@ -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, @@ -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 = [ @@ -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) @@ -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}'" @@ -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: @@ -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 @@ -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 @@ -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 = [i for i in range(nparts)] for name in sim.model_names: gwm = sim.get_model(name) if "concentration()" in gwm.output.methods(): @@ -1234,6 +1295,121 @@ def build_gwt_model(sim, gwtname, rch_package): diff = np.nansum(diff) if diff > 10.25: raise AssertionError( - "Difference between output arrays: " - f"{diff :.2f} greater than tolerance" + f"Difference between output arrays: {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") \ No newline at end of file From 7a824c5bd7e4b244878d530e977e210f9da7cceb Mon Sep 17 00:00:00 2001 From: jlarsen Date: Fri, 20 Dec 2024 17:02:09 -0800 Subject: [PATCH 3/5] linting --- autotest/test_model_splitter.py | 19 ++++++------------- 1 file changed, 6 insertions(+), 13 deletions(-) diff --git a/autotest/test_model_splitter.py b/autotest/test_model_splitter.py index 590c6f894..faff38bb1 100644 --- a/autotest/test_model_splitter.py +++ b/autotest/test_model_splitter.py @@ -1301,7 +1301,6 @@ def build_gwt_model(sim, gwtname, rch_package): @requires_exe("mf6") def test_timeseries(function_tmpdir): - sim = MFSimulation( sim_name="np001", sim_ws=function_tmpdir, @@ -1326,16 +1325,11 @@ def test_timeseries(function_tmpdir): preconditioner_drop_tolerance=0.01, number_orthogonalizations=2, ) - gwf = flopy.mf6.ModflowGwf(sim,) + gwf = flopy.mf6.ModflowGwf( + sim, + ) dis = flopy.mf6.ModflowGwfdis( - gwf, - nlay=1, - nrow=1, - ncol=10, - delr=500, - delc=500, - top=100, - botm=50 + gwf, nlay=1, nrow=1, ncol=10, delr=500, delc=500, top=100, botm=50 ) ic = flopy.mf6.ModflowGwfic(gwf, strt=95) @@ -1377,7 +1371,7 @@ def test_timeseries(function_tmpdir): "filename": "drn_ts.ts", "timeseries": [(0.0, 60.0), (100000.0, 60.0)], "time_series_namerecord": "drn_1", - "interpolation_methodrecord": "linearend" + "interpolation_methodrecord": "linearend", } drn = flopy.mf6.ModflowGwfdrn( gwf, @@ -1402,7 +1396,6 @@ def test_timeseries(function_tmpdir): sim.write_simulation() sim.run_simulation() - mfs = Mf6Splitter(sim) mask = mfs.optimize_splitting_mask(2) new_sim = mfs.split_model(mask) @@ -1412,4 +1405,4 @@ def test_timeseries(function_tmpdir): success, _ = new_sim.run_simulation() if not success: - raise AssertionError("Timeseries split simulation did not properly run") \ No newline at end of file + raise AssertionError("Timeseries split simulation did not properly run") From 7e7d0be64a2763631210ec8ccdd3083a7ca68547 Mon Sep 17 00:00:00 2001 From: jlarsen Date: Fri, 20 Dec 2024 17:03:50 -0800 Subject: [PATCH 4/5] Linting 2 --- autotest/test_model_splitter.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/autotest/test_model_splitter.py b/autotest/test_model_splitter.py index faff38bb1..d0f1b9dee 100644 --- a/autotest/test_model_splitter.py +++ b/autotest/test_model_splitter.py @@ -1264,7 +1264,7 @@ def build_gwt_model(sim, gwtname, rch_package): new_sim.run_simulation() # compare results for each of the models - splits = [i for i in range(nparts)] + splits = list(i for i in range(nparts)) for name in sim.model_names: gwm = sim.get_model(name) if "concentration()" in gwm.output.methods(): @@ -1295,7 +1295,8 @@ def build_gwt_model(sim, gwtname, rch_package): diff = np.nansum(diff) if diff > 10.25: raise AssertionError( - f"Difference between output arrays: {diff :.2f} greater than tolerance" + f"Difference between output arrays: " + f"{diff :.2f} greater than tolerance" ) From 1bb508fba9976b60b535873033884f4e926dbd4c Mon Sep 17 00:00:00 2001 From: jlarsen Date: Mon, 23 Dec 2024 09:09:47 -0800 Subject: [PATCH 5/5] linting --- autotest/test_model_splitter.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/autotest/test_model_splitter.py b/autotest/test_model_splitter.py index d0f1b9dee..a05708ec4 100644 --- a/autotest/test_model_splitter.py +++ b/autotest/test_model_splitter.py @@ -1264,7 +1264,7 @@ def build_gwt_model(sim, gwtname, rch_package): new_sim.run_simulation() # compare results for each of the models - splits = list(i for i in range(nparts)) + splits = range(nparts) for name in sim.model_names: gwm = sim.get_model(name) if "concentration()" in gwm.output.methods():