From c877fe5ea2936d4139e8675ee7f1084ffd85632c Mon Sep 17 00:00:00 2001 From: Joshua Larsen Date: Mon, 14 Aug 2023 08:06:15 -0700 Subject: [PATCH] updates(Mf6Splitter): control record and additional splitting checks * added checks for hfb and inactive models * updated remap_array and remap_mflist to maintain constants, external ascii files, and binary external files --- autotest/test_model_splitter.py | 128 ++++++++++++++++++++ flopy/mf6/utils/model_splitter.py | 190 +++++++++++++++++++++++++++++- 2 files changed, 312 insertions(+), 6 deletions(-) diff --git a/autotest/test_model_splitter.py b/autotest/test_model_splitter.py index 08cbb42c5..621bbaab7 100644 --- a/autotest/test_model_splitter.py +++ b/autotest/test_model_splitter.py @@ -1,4 +1,5 @@ import numpy as np +import flopy import pytest from autotest.conftest import get_example_data_path from modflow_devtools.markers import requires_exe, requires_pkg @@ -245,3 +246,130 @@ def test_save_load_node_mapping(function_tmpdir): new_heads = mfsplit2.reconstruct_array(array_dict) err_msg = "Heads from original and split models do not match" np.testing.assert_allclose(new_heads, original_heads, err_msg=err_msg) + + +def test_control_records(function_tmpdir): + nrow = 10 + ncol = 10 + nper = 3 + + sim = flopy.mf6.MFSimulation() + ims = flopy.mf6.ModflowIms(sim, complexity="SIMPLE") + + tdis = flopy.mf6.ModflowTdis( + sim, + nper=nper, + perioddata=( + (1., 1, 1.), + (1., 1, 1.), + (1., 1, 1.) + ) + ) + + gwf = flopy.mf6.ModflowGwf(sim, save_flows=True) + + botm2 = np.ones((nrow, ncol)) * 20 + dis = flopy.mf6.ModflowGwfdis( + gwf, + nlay=2, + nrow=nrow, + ncol=ncol, + delr=1, + delc=1, + top=35, + botm=[30, botm2], + idomain=1 + ) + + ic = flopy.mf6.ModflowGwfic(gwf, strt=32) + npf = flopy.mf6.ModflowGwfnpf( + gwf, + k=[ + 1., + { + "data": np.ones((10, 10)) * 0.75, + "filename": "k.l2.txt", + "iprn": 1, + "factor": 1 + } + ], + k33=[ + np.ones((nrow, ncol)), + { + "data": np.ones((nrow, ncol)) * 0.5, + "filename": "k33.l2.bin", + "iprn": 1, + "factor": 1, + "binary": True + } + ] + + ) + + wel_rec = [((0, 4, 5), -10), ] + + spd = { + 0: wel_rec, + 1: { + "data": wel_rec, + "filename": "wel.1.txt" + }, + 2: { + "data": wel_rec, + "filename": "wel.2.bin", + "binary": True + } + } + + wel = flopy.mf6.ModflowGwfwel( + gwf, + stress_period_data=spd + ) + + chd_rec = [] + for cond, j in ((30, 0), (22, 9)): + for i in range(10): + chd_rec.append(((0, i, j), cond)) + + chd = flopy.mf6.ModflowGwfchd( + gwf, + stress_period_data={0: chd_rec} + ) + + arr = np.zeros((10, 10), dtype=int) + arr[0:5, :] = 1 + + mfsplit = flopy.mf6.utils.Mf6Splitter(sim) + new_sim = mfsplit.split_model(arr) + + ml1 = new_sim.get_model("model_1") + + kls = ml1.npf.k._data_storage.layer_storage.multi_dim_list + if kls[0].data_storage_type.value != 2: + raise AssertionError("Constants not being preserved for MFArray") + + if kls[1].data_storage_type.value != 3 or kls[1].binary: + raise AssertionError( + "External ascii files not being preserved for MFArray" + ) + + k33ls = ml1.npf.k33._data_storage.layer_storage.multi_dim_list + if k33ls[1].data_storage_type.value != 3 or not k33ls[1].binary: + raise AssertionError( + "Binary file input not being preserved for MFArray" + ) + + spd_ls1 = ml1.wel.stress_period_data._data_storage[1].layer_storage.multi_dim_list[0] + spd_ls2 = ml1.wel.stress_period_data._data_storage[2].layer_storage.multi_dim_list[0] + + if spd_ls1.data_storage_type.value != 3 or spd_ls1.binary: + raise AssertionError( + "External ascii files not being preserved for MFList" + ) + + if spd_ls2.data_storage_type.value != 3 or not spd_ls2.binary: + raise AssertionError( + "External binary file input not being preseved for MFList" + ) + + diff --git a/flopy/mf6/utils/model_splitter.py b/flopy/mf6/utils/model_splitter.py index 6b1905dec..90ded666f 100644 --- a/flopy/mf6/utils/model_splitter.py +++ b/flopy/mf6/utils/model_splitter.py @@ -1,5 +1,6 @@ import inspect import json +import warnings import numpy as np @@ -366,6 +367,7 @@ def optimize_splitting_mask(self, nparts): adv_pkg_weights = np.zeros((ncpl,), dtype=int) lak_array = np.zeros((ncpl,), dtype=int) laks = [] + hfbs = [] for _, package in self._model.package_dict.items(): if isinstance( package, @@ -373,8 +375,13 @@ def optimize_splitting_mask(self, nparts): modflow.ModflowGwfsfr, modflow.ModflowGwfuzf, modflow.ModflowGwflak, + modflow.ModflowGwfhfb ), ): + if isinstance(package, modflow.ModflowGwfhfb): + hfbs.append(package) + continue + if isinstance(package, modflow.ModflowGwflak): cellids = package.connectiondata.array.cellid else: @@ -394,6 +401,8 @@ def optimize_splitting_mask(self, nparts): else: adv_pkg_weights[nodes] += 1 + + for nn, neighbors in self._modelgrid.neighbors().items(): weight = np.count_nonzero(idomain[:, nn]) adv_weight = adv_pkg_weights[nn] @@ -410,6 +419,23 @@ def optimize_splitting_mask(self, nparts): mnum = np.unique(membership[idx])[0] membership[idx] = mnum + if hfbs: + for hfb in hfbs: + for recarray in hfb.stress_period_data.data.values(): + cellids1 = recarray.cellid1 + cellids2 = recarray.cellid2 + _, nodes1 = self._cellid_to_layer_node(cellids1) + _, nodes2 = self._cellid_to_layer_node(cellids2) + mnums1 = membership[nodes1] + mnums2 = membership[nodes2] + ev = np.equal(mnums1, mnums2) + if np.alltrue(ev): + continue + idx = np.where(~ev)[0] + mnum_to = mnums1[idx] + adj_nodes = nodes2[idx] + membership[adj_nodes] = mnum_to + return membership.reshape(shape) def reconstruct_array(self, arrays): @@ -616,6 +642,24 @@ def _remap_nodes(self, array): """ array = np.ravel(array) + idomain = self._modelgrid.idomain.reshape((-1, self._ncpl)) + mkeys = np.unique(array) + bad_keys = [] + for mkey in mkeys: + count = 0 + mask = array[array == mkey] + for arr in idomain: + check = arr[mask] + count += np.count_nonzero(check) + if count == 0: + bad_keys.append(mkey) + + if bad_keys: + raise KeyError( + f"{bad_keys} are not in the active model extent; " + f"please adjust the model splitting array" + ) + if self._modelgrid.iverts is None: self._map_iac_ja_connections() else: @@ -1012,7 +1056,53 @@ def _remap_disu(self, mapped_data): return mapped_data - def _remap_array(self, item, mfarray, mapped_data): + def _remap_transient_array(self, item, mftransient, mapped_data): + """ + Method to split and remap transient arrays to new models + + Parameters + ---------- + item : str + variable name + mftransient : mfdataarray.MFTransientArray + transient array object + mapped_data : dict + dictionary of remapped package data + + Returns + ------- + dict + """ + if mftransient.array is None: + return 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 + ] + + 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, values in d0.items(): + mapped_data[mkey][item] = values + + return mapped_data + + def _remap_array(self, item, mfarray, mapped_data, **kwargs): """ Method to remap array nodes to each model @@ -1029,10 +1119,25 @@ def _remap_array(self, item, mfarray, mapped_data): ------- dict """ - if mfarray.array is None: - return mapped_data - + how = kwargs.pop("how", []) + binary = kwargs.pop("binary", []) + fnames = kwargs.pop("fnames", None) if not hasattr(mfarray, "size"): + if mfarray.array is None: + return mapped_data + + how = [ + i.data_storage_type.value for i in + mfarray._data_storage.layer_storage.multi_dim_list + ] + binary = [ + i.binary for i in + mfarray._data_storage.layer_storage.multi_dim_list + ] + fnames = [ + i.fname for i in + mfarray._data_storage.layer_storage.multi_dim_list + ] mfarray = mfarray.array nlay = 1 @@ -1068,11 +1173,44 @@ def _remap_array(self, item, mfarray, mapped_data): new_array[new_nodes] = original_arr[old_nodes] + if how: + new_input = [] + i0 = 0 + i1 = new_ncpl + lay = 0 + for h in how: + if h == 1: + # internal array + new_input.append(new_array[i0: i1]) + elif h == 2: + # constant, parse the original array data + new_input.append(original_arr[ncpl * lay]) + else: + # external array + tmp = fnames[lay].split(".") + filename = f"{'.'.join(tmp[:-1])}.{mkey}.{tmp[-1]}" + + cr = { + "filename": filename, + "factor": 1, + "iprn": 1, + "data": new_array[i0: i1], + "binary": binary[lay] + } + + new_input.append(cr) + + i0 += new_ncpl + i1 += new_ncpl + lay += 1 + + new_array = new_input + mapped_data[mkey][item] = new_array return mapped_data - def _remap_mflist(self, item, mflist, mapped_data, transient=False): + def _remap_mflist(self, item, mflist, mapped_data, transient=False, **kwargs): """ Method to remap mflist data to each model @@ -1087,15 +1225,21 @@ def _remap_mflist(self, item, mflist, mapped_data, transient=False): transient : bool flag to indicate this is transient stress period data flag is needed to trap for remapping mover data. + Returns ------- dict """ mvr_remap = {} + how = kwargs.pop("how", 1) + binary = kwargs.pop("binary", False) + fname = kwargs.pop("fname", None) if hasattr(mflist, "array"): if mflist.array is None: return mapped_data recarray = mflist.array + how = mflist._data_storage._data_storage_type.value + binary = mflist._data_storage.layer_storage.multi_dim_list[0].binary else: recarray = mflist @@ -1123,6 +1267,16 @@ def _remap_mflist(self, item, mflist, mapped_data, transient=False): new_recarray = recarray[idx] new_recarray["cellid"] = new_cellids + if how == 3 and new_recarray is not None: + tmp = fname.split(".") + filename = f"{'.'.join(tmp[:-1])}.{mkey}.{tmp[-1]}" + + new_recarray = { + "data": new_recarray, + "binary": binary, + "filename": filename + } + mapped_data[mkey][item] = new_recarray if not transient: @@ -2539,9 +2693,28 @@ def _remap_transient_list(self, item, mftransientlist, mapped_data): mapped_data[mkey][item] = None return mapped_data + how = { + p: i.layer_storage.multi_dim_list[0].data_storage_type.value for p, i in + mftransientlist._data_storage.items() + } + binary = { + p: i.layer_storage.multi_dim_list[0].binary for p, i in + mftransientlist._data_storage.items() + } + fnames = { + p: i.layer_storage.multi_dim_list[0].fname for p, i in + mftransientlist._data_storage.items() + } + for per, recarray in mftransientlist.data.items(): d, mvr_remaps = self._remap_mflist( - item, recarray, mapped_data, transient=True + item, + recarray, + mapped_data, + transient=True, + how=how[per], + binary=binary[per], + fname=fnames[per] ) for mkey in self._model_dict.keys(): if mapped_data[mkey][item] is None: @@ -2714,6 +2887,11 @@ def _remap_package(self, package, ismvr=False): for k, v in self._offsets[mkey].items(): mapped_data[mkey][k] = v + elif isinstance(value, mfdataarray.MFTransientArray): + mapped_data = self._remap_transient_array( + item, value, mapped_data + ) + elif isinstance(value, mfdataarray.MFArray): mapped_data = self._remap_array(item, value, mapped_data)