diff --git a/autotest/test_model_splitter.py b/autotest/test_model_splitter.py index fae6dbcb9..8b2e1c5d9 100644 --- a/autotest/test_model_splitter.py +++ b/autotest/test_model_splitter.py @@ -1,5 +1,6 @@ import numpy as np import pytest +import yaml from modflow_devtools.markers import requires_exe, requires_pkg from modflow_devtools.misc import set_dir @@ -570,14 +571,14 @@ def test_transient_array(function_tmpdir): 0, 2, ): - d[key] = g.sto.steady_state.get_data(key).get_data() + d[key] = g.sto.steady_state.get_data(key) 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() + d[key] = g.sto.transient.get_data(key) assert d == transient, ( "storage package transient dictionary " + f"does not match for model '{name}'" @@ -859,3 +860,467 @@ def test_unstructured_complex_disu(function_tmpdir): diff = np.abs(heads - new_heads) if np.max(diff) > 1e-07: raise AssertionError("Reconstructed head results outside of tolerance") + + +@requires_exe("mf6") +@requires_pkg("pymetis") +@requires_pkg("scipy") +def test_multi_model(function_tmpdir): + from scipy.spatial import KDTree + + def string2geom(geostring, conversion=None): + if conversion is None: + multiplier = 1.0 + else: + multiplier = float(conversion) + res = [] + for line in geostring.split("\n"): + if not any(line): + continue + line = line.strip() + line = line.split(" ") + x = float(line[0]) * multiplier + y = float(line[1]) * multiplier + res.append((x, y)) + return res + + sim_path = function_tmpdir + split_sim_path = sim_path / "model_split" + data_path = get_example_data_path() + + ascii_file = data_path / "geospatial/fine_topo.asc" + fine_topo = flopy.utils.Raster.load(ascii_file) + + with open(data_path / "groundwater2023/geometries.yml") as foo: + geometry = yaml.safe_load(foo) + + Lx = 180000 + Ly = 100000 + dx = 2500.0 + dy = 2500.0 + nrow = int(Ly / dy) + 1 + ncol = int(Lx / dx) + 1 + boundary = string2geom(geometry["boundary"]) + bp = np.array(boundary) + + stream_segs = ( + geometry["streamseg1"], + geometry["streamseg2"], + geometry["streamseg3"], + geometry["streamseg4"], + ) + sgs = [string2geom(sg) for sg in stream_segs] + + modelgrid = flopy.discretization.StructuredGrid( + nlay=1, + delr=np.full(ncol, dx), + delc=np.full(nrow, dy), + xoff=0.0, + yoff=0.0, + top=np.full((nrow, ncol), 1000.0), + botm=np.full((1, nrow, ncol), -100.0), + ) + + ixs = flopy.utils.GridIntersect(modelgrid, method="vertex", rtree=True) + result = ixs.intersect( + [ + boundary, + ], + shapetype="Polygon", + ) + r, c = list(zip(*list(result.cellids))) + idomain = np.zeros(modelgrid.shape, dtype=int) + idomain[:, r, c] = 1 + modelgrid._idomain = idomain + + top = fine_topo.resample_to_grid( + modelgrid, + band=fine_topo.bands[0], + method="linear", + extrapolate_edges=True, + ) + modelgrid._top = top + + # intersect stream segments + cellids = [] + lengths = [] + for sg in stream_segs: + sg = string2geom(sg) + v = ixs.intersect(sg, shapetype="LineString", sort_by_cellid=True) + cellids += v["cellids"].tolist() + lengths += v["lengths"].tolist() + + r, c = list(zip(*cellids)) + idomain[:, r, c] = 2 + modelgrid._idomain = idomain + + nlay = 5 + dv0 = 5.0 + hyd_cond = 10.0 + hk = np.full((nlay, nrow, ncol), hyd_cond) + hk[1, :, 25:] = hyd_cond * 0.001 + hk[3, :, 10:] = hyd_cond * 0.00005 + + # drain leakage + leakance = hyd_cond / (0.5 * dv0) + + drn_data = [] + for cellid, length in zip(cellids, lengths): + x = modelgrid.xcellcenters[cellid] + width = 5.0 + (14.0 / Lx) * (Lx - x) + conductance = leakance * length * width + if not isinstance(cellid, tuple): + cellid = (cellid,) + drn_data.append((0, *cellid, top[cellid], conductance)) + + discharge_data = [] + area = dx * dy + for r in range(nrow): + for c in range(ncol): + if idomain[0, r, c] == 1: + conductance = leakance * area + discharge_data.append( + (0, r, c, top[r, c] - 0.5, conductance, 1.0) + ) + + topc = np.zeros((nlay, nrow, ncol), dtype=float) + botm = np.zeros((nlay, nrow, ncol), dtype=float) + topc[0] = modelgrid.top.copy() + botm[0] = topc[0] - dv0 + for idx in range(1, nlay): + dv0 *= 1.5 + topc[idx] = botm[idx - 1] + botm[idx] = topc[idx] - dv0 + + strt = np.tile([modelgrid.top], (nlay, 1, 1)) + idomain = np.tile( + [ + modelgrid.idomain[0], + ], + (5, 1, 1), + ) + + # setup recharge + dist_from_riv = 10000.0 + + grid_xx = modelgrid.xcellcenters + grid_yy = modelgrid.ycellcenters + riv_idxs = np.array(cellids) + riv_xx = grid_xx[riv_idxs[:, 0], riv_idxs[:, 1]] + riv_yy = grid_yy[riv_idxs[:, 0], riv_idxs[:, 1]] + + river_xy = np.column_stack((riv_xx.ravel(), riv_yy.ravel())) + grid_xy = np.column_stack((grid_xx.ravel(), grid_yy.ravel())) + tree = KDTree(river_xy) + distance, index = tree.query(grid_xy) + + index2d = index.reshape(nrow, ncol) + distance2d = distance.reshape(nrow, ncol) + + mountain_array = np.asarray(distance2d > dist_from_riv).nonzero() + mountain_idxs = np.array(list(zip(mountain_array[0], mountain_array[1]))) + + valley_array = np.asarray(distance2d <= dist_from_riv).nonzero() + valley_idxs = np.array(list(zip(valley_array[0], valley_array[1]))) + + max_recharge = 0.0001 + + rch_orig = max_recharge * np.ones((nrow, ncol)) + + rch_mnt = np.zeros((nrow, ncol)) + for idx in mountain_idxs: + rch_mnt[idx[0], idx[1]] = max_recharge + + rch_val = np.zeros((nrow, ncol)) + for idx in valley_idxs: + rch_val[idx[0], idx[1]] = max_recharge + + sim = flopy.mf6.MFSimulation( + sim_ws=sim_path, + exe_name="mf6", + memory_print_option="summary", + ) + + nper = 10 + nsteps = 1 + year = 365.25 + dt = 1000 * year + tdis = flopy.mf6.ModflowTdis( + sim, nper=nper, perioddata=nper * [(nsteps * dt, nsteps, 1.0)] + ) + + gwfname = "gwf" + + imsgwf = flopy.mf6.ModflowIms( + sim, + complexity="simple", + print_option="SUMMARY", + linear_acceleration="bicgstab", + outer_maximum=1000, + inner_maximum=100, + outer_dvclose=1e-4, + inner_dvclose=1e-5, + preconditioner_levels=2, + relaxation_factor=0.0, + filename=f"{gwfname}.ims", + ) + + gwf = flopy.mf6.ModflowGwf( + sim, + modelname=gwfname, + print_input=False, + save_flows=True, + newtonoptions="NEWTON UNDER_RELAXATION", + ) + + dis = flopy.mf6.ModflowGwfdis( + gwf, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=dx, + delc=dy, + idomain=idomain, + top=modelgrid.top, + botm=botm, + xorigin=0.0, + yorigin=0.0, + ) + + ic = flopy.mf6.ModflowGwfic(gwf, strt=strt) + + npf = flopy.mf6.ModflowGwfnpf( + gwf, + save_specific_discharge=True, + icelltype=1, + k=hk, + ) + + sto = flopy.mf6.ModflowGwfsto( + gwf, + save_flows=True, + iconvert=1, + ss=0.00001, + sy=0.35, + steady_state={0: True, 1: False}, + transient={0: False, 1: True}, + ) + + rch0 = flopy.mf6.ModflowGwfrcha( + gwf, + pname="rch_original", + recharge={0: rch_orig, 1: 0.0}, + filename="gwf_original.rch", + ) + + rch1 = flopy.mf6.ModflowGwfrcha( + gwf, + pname="rch_mountain", + recharge={1: rch_mnt}, + auxiliary="CONCENTRATION", + aux={1: 1.0}, + filename="gwf_mountain.rch", + ) + + rch2 = flopy.mf6.ModflowGwfrcha( + gwf, + pname="rch_valley", + recharge={1: rch_val}, + auxiliary="CONCENTRATION", + aux={1: 1.0}, + filename="gwf_valley.rch", + ) + + drn = flopy.mf6.ModflowGwfdrn( + gwf, + stress_period_data=drn_data, + pname="river", + filename=f"{gwfname}_riv.drn", + ) + + drn_gwd = flopy.mf6.ModflowGwfdrn( + gwf, + auxiliary=["depth"], + auxdepthname="depth", + stress_period_data=discharge_data, + pname="gwd", + filename=f"{gwfname}_gwd.drn", + ) + + wel_spd = {0: [[4, 20, 30, 0.0], [2, 20, 60, 0.0], [2, 30, 50, 0.0]]} + + wel = flopy.mf6.ModflowGwfwel( + gwf, + print_input=False, + print_flows=False, + stress_period_data=wel_spd, + ) + + oc = flopy.mf6.ModflowGwfoc( + gwf, + head_filerecord=f"{gwf.name}.hds", + budget_filerecord=f"{gwf.name}.cbc", + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + printrecord=[("BUDGET", "ALL")], + ) + + sim.register_ims_package(imsgwf, [gwf.name]) + + def build_gwt_model(sim, gwtname, rch_package): + conc_start = 0.0 + diffc = 0.0 + alphal = 0.1 + porosity = 0.35 + gwf = sim.get_model("gwf") + modelgrid = gwf.modelgrid + + gwt = flopy.mf6.ModflowGwt( + sim, + modelname=gwtname, + print_input=False, + save_flows=True, + ) + + nlay, nrow, ncol = modelgrid.shape + + dis = flopy.mf6.ModflowGwtdis( + gwt, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=dx, + delc=dy, + idomain=modelgrid.idomain, + top=modelgrid.top, + botm=botm, + xorigin=0.0, + yorigin=0.0, + ) + + # initial conditions + ic = flopy.mf6.ModflowGwtic( + gwt, strt=conc_start, filename=f"{gwtname}.ic" + ) + + # advection + adv = flopy.mf6.ModflowGwtadv( + gwt, scheme="tvd", filename=f"{gwtname}.adv" + ) + + # dispersion + dsp = flopy.mf6.ModflowGwtdsp( + gwt, + diffc=diffc, + alh=alphal, + alv=alphal, + ath1=0.0, + atv=0.0, + filename=f"{gwtname}.dsp", + ) + + # mass storage and transfer + mst = flopy.mf6.ModflowGwtmst( + gwt, porosity=porosity, filename=f"{gwtname}.mst" + ) + + # sources + sourcerecarray = [ + (rch_package, "AUX", "CONCENTRATION"), + ] + ssm = flopy.mf6.ModflowGwtssm( + gwt, sources=sourcerecarray, filename=f"{gwtname}.ssm" + ) + + # output control + oc = flopy.mf6.ModflowGwtoc( + gwt, + budget_filerecord=f"{gwtname}.cbc", + concentration_filerecord=f"{gwtname}.ucn", + saverecord=[("CONCENTRATION", "ALL"), ("BUDGET", "ALL")], + ) + + return gwt + + imsgwt = flopy.mf6.ModflowIms( + sim, + complexity="complex", + print_option="SUMMARY", + linear_acceleration="bicgstab", + outer_maximum=1000, + inner_maximum=100, + outer_dvclose=1e-4, + inner_dvclose=1e-5, + filename="gwt.ims", + ) + + gwt_mnt = build_gwt_model(sim, "gwt_mnt", "rch_mountain") + sim.register_ims_package(imsgwt, [gwt_mnt.name]) + + gwt_val = build_gwt_model(sim, "gwt_val", "rch_valley") + sim.register_ims_package(imsgwt, [gwt_val.name]) + + gwfgwt = flopy.mf6.ModflowGwfgwt( + sim, + exgtype="GWF6-GWT6", + exgmnamea=gwfname, + exgmnameb=gwt_mnt.name, + filename="gwfgwt_mnt.exg", + ) + + gwfgwt = flopy.mf6.ModflowGwfgwt( + sim, + exgtype="GWF6-GWT6", + exgmnamea=gwfname, + exgmnameb=gwt_val.name, + filename="gwfgwt_val.exg", + ) + + sim.write_simulation() + sim.run_simulation() + + nparts = 2 + mfs = Mf6Splitter(sim) + array = mfs.optimize_splitting_mask(nparts) + new_sim = mfs.split_multi_model(array) + new_sim.set_sim_path(split_sim_path) + new_sim.write_simulation() + new_sim.run_simulation() + + # compare results for each of the models + splits = [i for i in range(nparts)] + for name in sim.model_names: + gwm = sim.get_model(name) + if "concentration()" in gwm.output.methods(): + X = gwm.output.concentration().get_alldata()[-1] + else: + X = gwm.output.head().get_alldata()[-1] + + array_dict = {} + for split in splits: + mname = f"{name}_{split}" + sp_gwm = new_sim.get_model(mname) + if "concentration()" in sp_gwm.output.methods(): + X0 = sp_gwm.output.concentration().get_alldata()[-1] + else: + X0 = sp_gwm.output.head().get_alldata()[-1] + + array_dict[split] = X0 + + X_split = mfs.reconstruct_array(array_dict) + + err_msg = ( + f"Outputs from {name} and split model " f"are not within tolerance" + ) + X_split[idomain == 0] = np.nan + X[idomain == 0] = np.nan + if name == "gwf": + np.testing.assert_allclose( + X, X_split, equal_nan=True, err_msg=err_msg + ) + else: + diff = np.abs(X_split - X) + diff = np.nansum(diff) + if diff > 10.25: + raise AssertionError( + f"Difference between output arrays: {diff :.2f} greater than tolerance" + ) diff --git a/flopy/mf6/utils/model_splitter.py b/flopy/mf6/utils/model_splitter.py index b678dce0b..e73a5cede 100644 --- a/flopy/mf6/utils/model_splitter.py +++ b/flopy/mf6/utils/model_splitter.py @@ -173,6 +173,9 @@ def __init__(self, sim, modelname=None): self._fdigits = 1 + # multi-model splitting attr + self._multimodel_exchange_gwf_names = {} + @property def new_simulation(self): """ @@ -941,6 +944,8 @@ def _create_sln_tdis(self): new_sim : MFSimulation object """ for pak in self._sim.sim_package_list: + if pak.package_abbr in ("gwfgwt", "gwfgwf", "gwfgwe"): + continue pak_cls = PackageContainer.package_factory(pak.package_abbr, "") signature = inspect.signature(pak_cls) d = {"simulation": self._new_sim, "loading_package": False} @@ -1019,6 +1024,7 @@ def _remap_filerecords(self, item, value, mapped_data): "budgetcsv_filerecord", "stage_filerecord", "obs_filerecord", + "concentration_filerecord", ): value = value.array if value is None: @@ -3009,7 +3015,13 @@ def _remap_package(self, package, ismvr=False): elif isinstance(value, mfdatascalar.MFScalarTransient): for mkey in self._model_dict.keys(): - mapped_data[mkey][item] = value._data_storage + val_dict = {} + for ( + perkey, + data_storage, + ) in value._data_storage.items(): + val_dict[perkey] = data_storage.get_data() + mapped_data[mkey][item] = val_dict elif isinstance(value, mfdatascalar.MFScalar): for mkey in self._model_dict.keys(): mapped_data[mkey][item] = value.data @@ -3072,28 +3084,29 @@ def _create_exchanges(self): dict """ d = {} + exchange_kwargs = {} built = [] nmodels = list(self._model_dict.keys()) - if self._model.name_file.newtonoptions is not None: - newton = self._model.name_file.newtonoptions.array - if isinstance(newton, list): - newton = True - else: - newton = None - - if self._model.npf.xt3doptions is not None: - xt3d = self._model.npf.xt3doptions.array - if isinstance(xt3d, list): - xt3d = True - else: - xt3d = None + if hasattr(self._model.name_file, "newtonoptions"): + if self._model.name_file.newtonoptions is not None: + newton = self._model.name_file.newtonoptions.array + if isinstance(newton, list): + exchange_kwargs["newton"] = True + + if hasattr(self._model, "npf"): + if self._model.npf.xt3doptions is not None: + xt3d = self._model.npf.xt3doptions.array + if isinstance(xt3d, list): + exchange_kwargs["xt3d"] = True if self._model_type.lower() == "gwf": extension = "gwfgwf" exchgcls = modflow.ModflowGwfgwf + check_multi_model = False elif self._model_type.lower() == "gwt": extension = "gwtgwt" exchgcls = modflow.ModflowGwtgwt + check_multi_model = True else: raise NotImplementedError() @@ -3111,6 +3124,14 @@ def _create_exchanges(self): if m1 == m0: continue exchange_data = [] + if check_multi_model: + if self._multimodel_exchange_gwf_names: + exchange_kwargs["gwfmodelname1"] = ( + self._multimodel_exchange_gwf_names[m0] + ) + exchange_kwargs["gwfmodelname2"] = ( + self._multimodel_exchange_gwf_names[m1] + ) for node0, exg_list in exg_nodes.items(): if grid0.idomain[node0] < 1: continue @@ -3145,9 +3166,9 @@ def _create_exchanges(self): exgmnameb=mname1, nexg=len(exchange_data), exchangedata=exchange_data, - filename=f"sim_{m0 :0{self._fdigits}d}_{m1 :0{self._fdigits}d}.{extension}", - newton=newton, - xt3d=xt3d, + filename=f"sim_{mname0}_{mname1}.{extension}", + pname=f"{mname0}_{mname1}", + **exchange_kwargs, ) d[f"{mname0}_{mname1}"] = exchg @@ -3155,7 +3176,8 @@ def _create_exchanges(self): for _, model in self._model_dict.items(): # turn off save_specific_discharge if it's on - model.npf.save_specific_discharge = None + if hasattr(model, "npf"): + model.npf.save_specific_discharge = None else: xc = self._modelgrid.xcellcenters.ravel() @@ -3164,17 +3186,27 @@ def _create_exchanges(self): for m0, model in self._model_dict.items(): exg_nodes = self._new_connections[m0]["external"] for m1 in nmodels: + exchange_data = [] if m1 in built: continue if m1 == m0: continue + + if check_multi_model: + if self._multimodel_exchange_gwf_names: + exchange_kwargs["gwfmodelname1"] = ( + self._multimodel_exchange_gwf_names[m0] + ) + exchange_kwargs["gwfmodelname2"] = ( + self._multimodel_exchange_gwf_names[m1] + ) + modelgrid0 = model.modelgrid modelgrid1 = self._model_dict[m1].modelgrid ncpl0 = modelgrid0.ncpl ncpl1 = modelgrid1.ncpl idomain0 = modelgrid0.idomain idomain1 = modelgrid1.idomain - exchange_data = [] for node0, exg_list in exg_nodes.items(): for exg in exg_list: if exg[0] != m1: @@ -3283,9 +3315,9 @@ def _create_exchanges(self): auxiliary=["ANGLDEGX", "CDIST"], nexg=len(exchange_data), exchangedata=exchange_data, - filename=f"sim_{m0 :0{self._fdigits}d}_{m1 :0{self._fdigits}d}.{extension}", - newton=newton, - xt3d=xt3d, + filename=f"sim_{mname0}_{mname1}.{extension}", + pname=f"{mname0}_{mname1}", + **exchange_kwargs, ) d[f"{mname0}_{mname1}"] = exchg @@ -3300,12 +3332,55 @@ def _create_exchanges(self): filename=f"{mname0}_{mname1}.mvr", ) - d[f"{mname0}_{mname1}_mvr"] = exchg + d[f"{mname0}_{mname1}_mvr"] = mvr built.append(m0) return d + def create_multi_model_exchanges(self, mname0, mname1): + """ + Method to create multi-model exchange packages, e.g., GWF-GWT + + Parameters + ---------- + mname0 : + mname1 : + + Returns + ------- + + """ + exchange_classes = { + "gwfgwt": modflow.ModflowGwfgwt, + "gwfgwe": modflow.ModflowGwfgwe, + } + ml0 = self._new_sim.get_model(mname0) + ml1 = self._new_sim.get_model(mname1) + + mtype0 = ml0.model_type[:3] + mtype1 = ml1.model_type[:3] + + if mtype0.lower() != "gwf": + raise AssertionError( + f"GWF must be the first specified type for multimodel " + f"exchanges: type supplied {mtype1.upper()} " + ) + + if mtype1.lower() not in ("gwt", "gwe"): + raise NotImplementedError( + f"Unsupported exchange type GWF-{mtype1.upper()}" + ) + + exchangecls = exchange_classes[f"{mtype0}{mtype1}"] + filename = f"{mname0}_{mname1}.exg" + exchangecls( + self._new_sim, + exgmnamea=mname0, + exgmnameb=mname1, + filename=filename, + ) + def split_model(self, array): """ User method to split a model based on an array @@ -3366,3 +3441,119 @@ def split_model(self, array): epaks = self._create_exchanges() return self._new_sim + + def split_multi_model(self, array): + """ + Method to split integrated models such as GWF-GWT or GWF-GWE models. + Note: this method will not work to split multiple connected GWF models + + Parameters + ---------- + array : np.ndarray + integer array of new model numbers. Array must either be of + dimension (NROW, NCOL), (NCPL), or (NNODES for unstructured grid + models). + + Returns + ------- + MFSimulation object + """ + if not self._allow_splitting: + raise AssertionError( + "Mf6Splitter cannot split a model that " + "is part of a split simulation" + ) + + # set number formatting string for file paths + array = np.array(array).astype(int) + s = str(np.max(array)) + self._fdigits = len(s) + + # get model names and types, assert that first model is a GWF model + model_names = self._sim.model_names + models = [self._sim.get_model(mn) for mn in model_names] + model_types = [type(ml) for ml in models] + + ix = model_types.index(modflow.ModflowGwf) + if ix != 0: + idxs = [ + ix, + ] + [idx for idx in range(len(model_names)) if idx != ix] + model_names = [model_names[idx] for idx in idxs] + models = [models[idx] for idx in idxs] + model_types = [model_types[idx] for idx in idxs] + self.switch_models(modelname=model_names[0], remap_nodes=True) + + # assert consistent idomain and modelgrid shapes! + shapes = [ml.modelgrid.shape for ml in models] + idomains = [ml.modelgrid.idomain for ml in models] + gwf_shape = shapes.pop(0) + gwf_idomain = idomains.pop(0) + + for ix, shape in enumerate(shapes): + idomain = idomains[ix] + mname = model_names[ix + 1] + if shape != gwf_shape: + raise AssertionError( + f"Model {mname} shape {shape} is not consistent with GWF " + f"model shape {gwf_shape}" + ) + + gwf_inactive = np.where(gwf_idomain == 0) + inactive = np.where(idomain == 0) + if not np.allclose(inactive, gwf_inactive): + raise AssertionError( + f"Model {mname} idomain is not consistent with GWF " + f"model idomain" + ) + + gwf_base = model_names[0] + model_labels = [ + f"{i :{self._fdigits}d}" for i in sorted(np.unique(array)) + ] + + self._multimodel_exchange_gwf_names = { + int(i): f"{gwf_base}_{i}" for i in model_labels + } + + new_sim = self.split_model(array) + for mname in model_names[1:]: + self.switch_models(modelname=mname, remap_nodes=False) + new_sim = self.split_model(array) + + for mbase in model_names[1:]: + for label in model_labels: + mname0 = f"{gwf_base}_{label}" + mname1 = f"{mbase}_{label}" + self.create_multi_model_exchanges(mname0, mname1) + + # register models to correct IMS package + solution_recarray = self._sim.name_file.solutiongroup.data[0] + sln_mname_cols = [ + i for i in solution_recarray.dtype.names if "slnmnames" in i + ] + if len(solution_recarray) > 1: + # need to associate solutions with solution groups + imspkgs = [] + imspkg_names = [] + for pkg in new_sim.sim_package_list: + if isinstance(pkg, modflow.ModflowIms): + imspkgs.append(pkg) + imspkg_names.append(pkg.filename) + + for record in solution_recarray: + fname = record.slnfname + ims_ix = imspkg_names.index(fname) + ims_obj = imspkgs[ims_ix] + mnames = [] + for col in sln_mname_cols: + mbase = record[col] + if mbase is None: + continue + + for label in model_labels: + mnames.append(f"{mbase}_{label}") + + new_sim.register_ims_package(ims_obj, mnames) + + return self._new_sim