diff --git a/docs/source/concepts/lgr.rst b/docs/source/concepts/lgr.rst index 13d81f57..bee938cc 100644 --- a/docs/source/concepts/lgr.rst +++ b/docs/source/concepts/lgr.rst @@ -2,18 +2,28 @@ Local grid refinement =========================================================== -In MODFLOW 6, two groundwater models can be tightly coupled in the same simulation, which allows for efficient "local grid refinement" (LGR; Mehl and others, 2006) and "semistructured" (Feinstein and others, 2016) configurations that combine multiple structured model layers at different resolutions (Langevin and others, 2017). Locally refined area(s) are conceptualized as separate model(s) that are linked to the surrounding (usually coarser) "parent" model through the GWF Model Exchange Package. Similarly, "semistructured" configurations are represented by multiple linked models for each layer or group of layers with the same resolution. +In MODFLOW 6, two groundwater models can be tightly coupled in the same simulation, which allows for efficient "local grid refinement" (LGR; Mehl and others, 2006) and "semistructured" (Feinstein and others, 2016) configurations that combine multiple structured model layers at different resolutions (Langevin and others, 2017). Locally refined areas are conceptualized as separate "child" models that are linked to the surrounding (usually coarser) "parent" model through the GWF Model Exchange (GWFGWF) Package. Similarly, "semistructured" configurations are represented by multiple linked models for each layer or group of layers with the same resolution. -Modflow-setup supports this functionality via an ``lgr:`` subblock within the ``setup_grid:`` block of the configuration file. The ``lgr:`` subblock consists of one of more subblocks, each keyed by a linked model name and containing configuration input for that linked model. Vertical layer refinement relative to the "parent" model, and the range of parent model layers that will be occupied by the linked "inset" model are also specified. +Modflow-setup supports local grid refinment via an ``lgr:`` subblock within the ``setup_grid:`` block of the configuration file. The ``lgr:`` subblock consists of one of more subblocks, each keyed by a linked model name and containing configuration input for that linked model. Vertical layer refinement relative to the "parent" model is also specified for each layer of the parent model. -For example, the following "parent" configuration for the :ref:`Pleasant Lake Example ` creates a local refinement model named "``pleasant_lgr_inset``" that spans layers 0 through 3 (the first 4 layers) of the parent model, at the same vertical resolution (1 inset model layer per parent model layer). The horizontal location and resolution of the inset model are specified in the inset model configuration file, in the same way that they are specified for any model. +For example, the following "parent" configuration for the :ref:`Pleasant Lake Example ` creates a local refinement model named "``pleasant_lgr_inset``" that spans all layers of the parent model, at the same vertical resolution (1 inset model layer per parent model layer). .. literalinclude:: ../../../examples/pleasant_lgr_parent.yml :language: yaml :start-at: lgr: :end-before: # Structured Discretization Package -Input from the ``lgr:`` subblock and the inset model configuration file(s) is passed to the :py:class:`Flopy Lgr Utility `, which helps create input for the GWF Model Exchange Package. +The horizontal location and resolution of the inset model are specified in the :ref:`inset model configuration file `, in the same way that they are specified for any model. In this example, the parent model has a uniform horizontal resolution of 200 meters, and the inset a uniform resolution of 40 meters (a horizontal refinement of 5 inset model cells per parent model cells). The inset model resolution must be a factor of the parent model resolution. + +.. image:: ../_static/pleasant_lgr.png + :width: 1200 + :alt: Example of LGR refinement in all layers. + +.. image:: ../_static/pleasant_lgr_xsection.png + :width: 1200 + :alt: Example of LGR refinement in all layers. + +Input from the ``lgr:`` subblock and the inset model configuration file(s) is passed to the :py:class:`Flopy Lgr Utility `, which helps create input for the GWF Model Exchange Package. Within the context of a Python session, inset model information is stored in a dictionary under an ``inset`` attribute attached to the parent model. For example, to access a Flopy model object for the above inset model from a parent model named ``model``: @@ -22,10 +32,58 @@ Within the context of a Python session, inset model information is stored in a d inset_model = model.inset['pleasant_lgr_inset'] + + +Specification of vertical refinement +----------------------------------------- +Vertical refinement in the LGR child grid is specified in the ``layer_refinement:`` item, as the number of child layers in each parent model layer. Currently vertical refinement is limited to even subdivision of parent model layers. Vertical refinement can be specified as an integer for uniform refinement across all parent model layers: + +.. code-block:: yaml + + layer_refinement: 1 + +a list with an entry for each parent layer: + +.. code-block:: yaml + + layer_refinement: [1, 1, 1, 0, 0] + +or a dictionary with entries for each parent layer that is refined: + +.. code-block:: yaml + + layer_refinement: + 0:1 + 1:1 + 2:1 + +Parent model layers with 0 specified refinement are excluded from the child model. The list and dictionary inputs above are equivalent, as unlisted layers in the dictionary are assigned default refinement values of 0. Refinement values > 1 result in even subdivision of the parent layers. Similar to one-way coupled inset models, LGR child model layers surfaces can be discretized at the finer child resolution from the original source data. In the example below, a 9-layer child model is set within the top 4 layers of a 5 layer parent model. The parent model ``lgr:`` configuration block is specified as: + +.. literalinclude:: ../../../mfsetup/tests/data/pleasant_vertical_lgr_parent.yml + :language: yaml + :start-at: lgr: + :end-before: # Structured Discretization Package + +In the child model ``dis:`` configuration block, raster surfaces that were used to define the parent model layer bottoms are specified at their respective desired locations within the child model grid; Modflow-setup then subdivides these to create the desired layer configuration. The layer specification in the child model ``dis:`` block must be consistent with the `layer_refinement:` input in the parent model configuration (see below). The image below shows + +.. literalinclude:: ../../../mfsetup/tests/data/pleasant_vertical_lgr_inset.yml + :language: yaml + :start-at: dis: + :end-before: # Recharge and Well packages are inherited + +The figure below shows a cross section through the model grid resulting from this configuration: + +.. image:: ../_static/pleasant_vlgr_xsection.png + :width: 1500 + :alt: Example of partial vertical LGR refinement with layer subdivision. + + **A few notes about LGR functionality in Modflow-setup** * **Locally refined "inset" models must be aligned with the parent model grid**, which also means that their horizontal resolution must be a factor of the "parent" model resolution. Modflow-setup handles the alignment automatically by "snapping" inset model grids to the lower left corner of the parent cell containing the lower left corner of the inset model grid (the inset model origin in real-world coordinates). -* Similarly, inset models need to align vertically with the parent model layers. Parent layers can be subdivided using the ``layer_refinement:`` input option. -* In principal, a **semistructured** configuration consisting of a "parent" model at the coarsest horizontal resolution, and one or more inset models representing [horizontally and potentially vertically] refined layers within the "parent" model domain is possible, but this configuration has not been tested with Modflow-setup. Please contact the Modflow-setup development team by posting a `GitHub issue `_ if you are interested in this. -* Similarly, multiple inset models at different horizontal locations, and even inset models within inset models should be achievable, but have not been tested. +* Similarly, inset models need to align vertically with the parent model layers. Parent layers can be subdivided using values > 1 in the ``layer_refinement:`` input option. +* Specifically, **the layer specification in the child model** ``dis:`` **block must be consistent with the** ``layer_refinement:`` **input in the parent model configuration**. For example, if a ``layer_refinement`` of ``3`` is specified for the last parent layer included in the child model domain, then the last two raster surfaces specified in the child model ``dis``: block must be specified with two layer bottoms in between. Similarly, the values in ``layer_refinement:`` in the parent model configuration must sum to ``nlay:`` specified in the child model ``dis:`` configuration block. +* Regardless of the supplied input, the child model bottom and parent model tops are aligned to remove any vertical gaps or overlap in the numerical model grid. If a raster surface is supplied for the child model bottom, the child bottom/parent top surface is based on the mean bottom elevations sampled to the child cells within each parent cell area. +* Child model ``layer_refinement:`` must start at the top of the parent model, and include a contiguous sequence of parent model layers. +* Multiple inset models at different horizontal locations, and even inset models within inset models should be achievable, but have not been tested. * **Multi-model configurations come with costs.** Each model within a MODFLOW 6 simulation carries its own set of files, including external text array and list input files to packages. As with a single model, when using the automated parameterization functionality in `pyEMU `_, the number of files is multiplied. At some point, at may be more efficient to work with individual models, and design the grids in such a way that boundary conditions along the model perimeters have a minimal impact on the predictions of interest. diff --git a/examples/Pleasant_lake_lgr_example.ipynb b/examples/Pleasant_lake_lgr_example.ipynb index 1a4b5760..8056d36c 100644 --- a/examples/Pleasant_lake_lgr_example.ipynb +++ b/examples/Pleasant_lake_lgr_example.ipynb @@ -277,6 +277,33 @@ "plt.colorbar(qmi)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Make a cross section of the grid" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(14, 5))\n", + "xs_line = [(553000, 390200), (558000, 390200)]\n", + "xs = flopy.plot.PlotCrossSection(model=m,\n", + " line={\"line\": xs_line}, ax=ax,\n", + " geographic_coords=True)\n", + "lc = xs.plot_grid(zorder=4)\n", + "xs2 = flopy.plot.PlotCrossSection(model=inset,\n", + " line={\"line\": xs_line}, ax=ax,\n", + " geographic_coords=True)\n", + "lc = xs2.plot_grid(zorder=4)\n", + "ax.set_ylim(190, 400)\n", + "plt.savefig('../../docs/source/_static/pleasant_lgr_xsection.png', bbox_inches='tight')" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -528,7 +555,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.2" + "version": "3.10.14" } }, "nbformat": 4, diff --git a/examples/pleasant_lgr_inset.yml b/examples/pleasant_lgr_inset.yml index 51ff9abc..ab163592 100644 --- a/examples/pleasant_lgr_inset.yml +++ b/examples/pleasant_lgr_inset.yml @@ -44,7 +44,7 @@ dis: length_units: 'meters' dimensions: # nrow and ncol are based on the buffer distance and grid spacing - nlay: 4 + nlay: 5 source_data: top: filename: 'data/pleasant/source_data/rasters/dem40m.tif' # DEM file; path relative to setup script diff --git a/examples/pleasant_lgr_parent.yml b/examples/pleasant_lgr_parent.yml index 70d8c719..a3df95b6 100644 --- a/examples/pleasant_lgr_parent.yml +++ b/examples/pleasant_lgr_parent.yml @@ -79,9 +79,12 @@ setup_grid: lgr: pleasant_lgr_inset: filename: 'pleasant_lgr_inset.yml' - layer_refinement: 1 # number of lgr model layers per parent model layer - parent_start_layer: 0 # uppermost parent model layer to include in lgr - parent_end_layer: 3 # lowermost parent model layer to include in lgr + # number of inset layers for each parent layer + # this can be an integer (vertically uniform refinement) + # a list (with an entry for each parent layer) + # or a dictionary (with entries for each parent layer that is refined; + # unlisted layers will default to zero) + layer_refinement: [1, 1, 1, 1, 1] # number of lgr model layers per parent model layer # Structured Discretization Package dis: diff --git a/mfsetup/bcs.py b/mfsetup/bcs.py index 680b7ec1..7869c2a8 100644 --- a/mfsetup/bcs.py +++ b/mfsetup/bcs.py @@ -394,7 +394,6 @@ def remove_inactive_bcs(pckg, external_files=False): Parameters ---------- - model : flopy model instance pckg : flopy package instance """ model = pckg.parent diff --git a/mfsetup/mf6model.py b/mfsetup/mf6model.py index 058ce76f..11edd56e 100644 --- a/mfsetup/mf6model.py +++ b/mfsetup/mf6model.py @@ -2,6 +2,7 @@ import os import shutil import time +from pathlib import Path import flopy import numpy as np @@ -36,7 +37,7 @@ ) from mfsetup.mfmodel import MFsetupMixin from mfsetup.mover import get_mover_sfr_package_input -from mfsetup.obs import setup_head_observations +from mfsetup.obs import remove_inactive_obs, setup_head_observations from mfsetup.oc import parse_oc_period_input from mfsetup.tdis import add_date_comments_to_tdis, setup_perioddata from mfsetup.units import convert_time_units @@ -200,6 +201,7 @@ def _set_idomain(self): if isinstance(self.lgr, dict): for k, v in self.lgr.items(): lgr_idomain[v.idomain == 0] = 0 + self._lgr_idomain2d = lgr_idomain[0] idomain_from_layer_elevations = make_idomain(self.dis.top.array, self.dis.botm.array, nodata=self._nodata_value, @@ -379,28 +381,55 @@ def create_lgr_models(self): parent_start_layer = v.get('parent_start_layer', 0) # parent_end_layer is specified as the last zero-based # parent layer that includes LGR refinement (not as a slice end) - parent_end_layer = v.get('parent_end_layer', self.nlay - 1) + 1 - ncppl = v.get('layer_refinement', 1) - specified_nlay_lgr = ncppl * (parent_end_layer - parent_start_layer) - specified_nlay_dis = inset_cfg['dis']['dimensions']['nlay'] - if specified_nlay_lgr != specified_nlay_dis: + parent_end_layer = v.get('parent_end_layer', self.nlay - 1) + # the layer refinement can be specified as an int, a list or a dict + ncppl_input = v.get('layer_refinement', 1) + if np.isscalar(ncppl_input): + ncppl = np.array([0] * self.modelgrid.nlay) + ncppl[parent_start_layer:parent_end_layer+1] = ncppl_input + elif isinstance(ncppl_input, list): + if not len(ncppl_input) == self.modelgrid.nlay: + raise ValueError( + "Configuration input: layer_refinement specified as" + "a list must include a value for every layer." + ) + ncppl = ncppl_input.copy() + elif isinstance(ncppl_input, dict): + ncppl = [ncppl_input.get(i, 0) for i in range(self.modelgrid.nlay)] + else: + raise ValueError("Configuration input: Unsupported input for " + "layer_refinement: supply an int, list or dict.") + + # refined layers must be consecutive, starting from layer 1 + is_refined = (np.array(ncppl) > 0).astype(int) + last_refined_layer = max(np.where(is_refined > 0)[0]) + consecutive = all(np.diff(is_refined)[:last_refined_layer] == 0) + if (is_refined[0] != 1) | (not consecutive): + raise ValueError("Configuration input: layer_refinement must" + "include consecutive sequence of layers, " + "starting with the top layer.") + # check the specified DIS package input is consistent + # with the specified layer_refinement + specified_nlay_dis = inset_cfg['dis']['dimensions'].get('nlay') + # skip this check if nlay hasn't been entered into the configuration file yet + if specified_nlay_dis and (np.sum(ncppl) != specified_nlay_dis): raise ValueError( - f"Parent model layers of {parent_start_layer} to {parent_end_layer} " - f"and layer refinement of {ncppl} implies {specified_nlay_lgr} inset model layers.\n" + f"Configuration input: layer_refinement of {ncppl} " + f"implies {is_refined.sum()} inset model layers.\n" f"{specified_nlay_dis} inset model layers specified in DIS package.") # mapping between parent and inset model layers # that is used for copying input from parent model - parent_inset_layer_mapping = dict() + inset_parent_layer_mapping = dict() inset_k = -1 - for parent_k in range(parent_start_layer, parent_end_layer): - for i in range(ncppl): + for parent_k, n_inset_lay in enumerate(ncppl): + for i in range(n_inset_lay): inset_k += 1 - parent_inset_layer_mapping[parent_k] = inset_k + inset_parent_layer_mapping[inset_k] = parent_k self.inset[inset_model.name].cfg['parent']['inset_layer_mapping'] =\ - parent_inset_layer_mapping + inset_parent_layer_mapping # create idomain indicating area of parent grid that is LGR lgr_idomain = make_lgr_idomain(self.modelgrid, self.inset[inset_model.name].modelgrid, - parent_start_layer, parent_end_layer) + ncppl) # inset model horizontal refinement from parent resolution refinement = self.modelgrid.delr[0] / self.inset[inset_model.name].modelgrid.delr[0] @@ -412,10 +441,55 @@ def create_lgr_models(self): self.dis.top.array, self.dis.botm.array, lgr_idomain, ncpp, ncppl) inset_model._perioddata = self.perioddata - self._set_idomain() + # set parent model top in LGR area to bottom of LGR area + # this is an initial draft; + # bottom elevations are readjusted in sourcedata.py + # when inset model DIS package botm array is set up + # (set to mean of inset model bottom elevations + # within each parent cell) + # number of layers in parent model with LGR + n_parent_lgr_layers = np.sum(np.array(ncppl) > 0) + lgr_area = self.lgr[inset_model.name].idomain == 0 + self.dis.top[lgr_area[0]] =\ + self.lgr[inset_model.name].botmp[n_parent_lgr_layers -1][lgr_area[0]] + # set parent model layers in LGR area to zero-thickness + new_parent_botm = self.dis.botm.array.copy() + for k in range(n_parent_lgr_layers): + new_parent_botm[k][lgr_area[0]] = self.dis.top[lgr_area[0]] + self.dis.botm = new_parent_botm + self._update_top_botm_external_files() + + + def _update_top_botm_external_files(self): + """Update the external files after assigning new elevations to the + Discretization Package top and botm arrays; adjust idomain as needed.""" + # reset the model top + # (this step may not be needed if the "original top" functionality + # is limited to cases where there is a lake package, + # or if the "original top"/"lake bathymetry" functionality is eliminated + # and we instead require the top to be pre-processed) + original_top_file = Path(self.external_path, + f"{self.name}_{self.cfg['dis']['top_filename_fmt']}.original") + original_top_file.unlink(missing_ok=True) + self._setup_array('dis', 'top', + data={0: self.dis.top.array}, + datatype='array2d', resample_method='linear', + write_fmt='%.2f', dtype=float) + # _set_idomain() regerates external files for bottom array + self._set_idomain() + def setup_lgr_exchanges(self): for inset_name, inset_model in self.inset.items(): + + # update cell information for computing any bottom exchanges + self.lgr[inset_name].top = inset_model.dis.top.array + self.lgr[inset_name].botm = inset_model.dis.botm.array + # update only the layers of the parent model below the child model + parent_top_below_child = np.sum(self.lgr[inset_name].ncppl > 0) -1 + self.lgr[inset_name].botmp[parent_top_below_child:] =\ + self.dis.botm.array[parent_top_below_child:] + # get the exchange data exchangelist = self.lgr[inset_name].get_exchange_data(angldegx=True, cdist=True) @@ -426,7 +500,9 @@ def setup_lgr_exchanges(self): # unpack the cellids and get their respective ibound values k1, i1, j1 = zip(*exchangedf['cellidm1']) k2, i2, j2 = zip(*exchangedf['cellidm2']) + # limit connections to active1 = self.idomain[k1, i1, j1] >= 1 + active2 = inset_model.idomain[k2, i2, j2] >= 1 # screen out connections involving an inactive cell @@ -454,6 +530,7 @@ def setup_lgr_exchanges(self): kwargs = get_input_arguments(kwargs, mf6.ModflowGwfgwf) gwfe = mf6.ModflowGwfgwf(self.simulation, **kwargs) + def setup_dis(self, **kwargs): """""" package = 'dis' @@ -624,7 +701,6 @@ def setup_rch(self, **kwargs): resample_method='nearest', write_fmt='%.6e', write_nodata=0.) - kwargs = self.cfg[package].copy() kwargs.update(self.cfg[package]['options']) kwargs = get_input_arguments(kwargs, mf6.ModflowGwfrcha) @@ -901,6 +977,9 @@ def write_input(self): external_files = self.cfg[pckg.lower()]['stress_period_data'] remove_inactive_bcs(package_instance, external_files=external_files) + if hasattr(self, 'obs'): + for obs_package_instance in self.obs: + remove_inactive_obs(obs_package_instance) # write the model with flopy # but skip the sfr package diff --git a/mfsetup/mfmodel.py b/mfsetup/mfmodel.py index 1b7e030c..3829c88d 100644 --- a/mfsetup/mfmodel.py +++ b/mfsetup/mfmodel.py @@ -120,6 +120,7 @@ def __init__(self, parent): self.inset = None # dictionary of inset models attached to LGR parent self._is_lgr = False # flag for lgr inset models self.lgr = None # holds flopy Lgr utility object + self._lgr_idomain2d = None # array of Lgr inset model locations within parent grid self.tmr = None # holds TMR class instance for TMR-type perimeter boundaries self._load = False # whether model is being made or loaded from existing files self.lake_info = None @@ -1301,12 +1302,13 @@ def _setup_basic_stress_package(self, package, flopy_package_class, if 'source_data' in kwargs: if package == 'wel': dropped_wells_file =\ - kwargs['output_files']['dropped_wells_file'].format(self.name) + kwargs.get('output_files', {})\ + .get('dropped_wells_file', '{}_dropped_wells.csv').format(self.name) df_sd = setup_wel_data(self, source_data=kwargs['source_data'], dropped_wells_file=dropped_wells_file) else: - df_sd = setup_basic_stress_data(self, **kwargs['source_data'], **kwargs['mfsetup_options']) + df_sd = setup_basic_stress_data(self, **kwargs['source_data'], **kwargs.get('mfsetup_options', dict())) if df_sd is not None and len(df_sd) > 0: dfs.append(df_sd) # set up package from parent model @@ -1514,7 +1516,13 @@ def setup_sfr(self, **kwargs): # create isfr array (where SFR cells will be populated) if self.version == 'mf6': active_cells = np.sum(self.idomain >= 1, axis=0) > 0 - #active_cells = self.idomain.sum(axis=0) > 0 + # For models with LGR, set the LGR area to isfr=0 + # to prevent SFR from being generated within the LGR area + # needed for LGR models that only have refinement + # in some layers (in other words, active parent model cells + # below the LGR inset) + if self.lgr: + active_cells[self._lgr_idomain2d == 0] = 0 else: active_cells = np.sum(self.ibound >= 1, axis=0) > 0 #active_cells = self.ibound.sum(axis=0) > 0 diff --git a/mfsetup/obs.py b/mfsetup/obs.py index 6868cf5b..605298e9 100644 --- a/mfsetup/obs.py +++ b/mfsetup/obs.py @@ -208,3 +208,28 @@ def make_obsname(name, unique_names={}, if name[slc] not in unique_names: return name[slc] return name[-maxlen:] + + +def remove_inactive_obs(obs_package_instance): + """Remove boundary conditions from cells that are inactive. + + Parameters + ---------- + obs_package_instance : flopy Modflow-6 Observation package instance + """ + model = obs_package_instance.parent + idomain = model.dis.idomain.array + if model.version != 'mf6': + raise NotImplementedError( + "obs.py::remove_inactive_obs(): " + "Support for removing MODFLOW 2005 not implemented.") + for obsfile, recarray in obs_package_instance.continuous.data.items(): + try: + k, i, j = zip(*recarray['id']) + except: + # for now, only support obs defined by k, i, j cellids + return + # cull any observations in inactive cells + is_active = idomain[k, i, j] > 0 + # update the flopy dataset + obs_package_instance.continuous.set_data({obsfile: recarray[is_active]}) diff --git a/mfsetup/sourcedata.py b/mfsetup/sourcedata.py index 05b268d1..09857608 100644 --- a/mfsetup/sourcedata.py +++ b/mfsetup/sourcedata.py @@ -502,7 +502,10 @@ def get_data(self): weight0=weight0) # repeat last layer if length of data is less than number of layers - if self.datatype == 'array3d' and i < (self.dest_model.nlay - 1): + # we probably don't want to do this for layer bottoms, + # but may want to for aquifer properties + if (self.variable != 'botm') and (self.datatype == 'array3d')\ + and (i < (self.dest_model.nlay - 1)): for j in range(i, self.dest_model.nlay): data[j] = data[i] @@ -1425,28 +1428,18 @@ def setup_array(model, package, var, data=None, model.cfg['dis']['top_filename_fmt']) top = model.load_array(model.cfg['dis']['griddata']['top'][0]['filename']) - # special case of LGR models - # with bottom connections to underlying parent cells - if model.version == 'mf6': - # (if model is an lgr inset model) - if model._is_lgr: - # regardless of what is specified for inset model bottom - # use top elevations of underlying parent model cells - nlay = model.cfg['dis']['dimensions']['nlay'] - if (nlay < model.parent.modelgrid.nlay): - # use the parent model bottoms - # mapped to the inset model grid by the Flopy Lgr util - lgr = model.parent.lgr[model.name] # Flopy Lgr inst. - data[nlay-1] = lgr.botm[nlay-1] - # set parent model top to top of - # first active parent model layer below this lgr domain - kp = lgr.get_parent_indices(nlay-1, 0, 0)[0] - model.parent.dis.top = lgr.botmp[kp] - j=2 - # fill missing layers if any if len(data) < model.nlay: all_surfaces = np.zeros((model.nlay + 1, model.nrow, model.ncol), dtype=float) * np.nan + # for LGR models, populate any missing bottom layer(s) + # from the refined parent grid + # this allows for the LGR model to be + # bounded on the bottom by a sub-divided parent layer surface + last_specified_layer = max(data.keys()) + if model._is_lgr and last_specified_layer < (model.nlay - 1): + #for k in range(last_specified_layer + 1, model.nlay): + # all_surfaces[k+1:] = model.parent.lgr[model.name].botm[k] + all_surfaces[model.nlay] = model.parent.lgr[model.name].botm[model.nlay -1] all_surfaces[0] = top for k, botm in data.items(): all_surfaces[k + 1] = botm @@ -1476,8 +1469,8 @@ def setup_array(model, package, var, data=None, # (so they won't get expanded again by fix_model_layer_conflicts) # only do this for mf6, where pinched out cells are allowed min_thickness = model.cfg['dis'].get('minimum_layer_thickness', 1) -# if model.version == 'mf6': -# botm[botm >= (top - min_thickness)] = np.nan + if model.version == 'mf6': + botm[botm >= (top - min_thickness)] = np.nan #for k, kbotm in enumerate(botm): # inlayer = lake_botm_elevations > kbotm[bathy != 0] @@ -1512,26 +1505,74 @@ def setup_array(model, package, var, data=None, else: model.dis.top = model.cfg['dis']['top'][0] data = {i: arr for i, arr in enumerate(botm)} - elif var in ['rech', 'recharge'] and simulate_high_k_lakes: - for per in range(model.nper): - if per == 0 and per not in data: - raise KeyError("No recharge input specified for first stress period.") - if per in data: - # assign high-k lake recharge for stress period - # only assign if precip and open water evaporation data were read - # (otherwise keep original values in recharge array) - last_data_array = data[per].copy() - if model.high_k_lake_recharge is not None: - data[per][model.isbc[0] == 2] = model.high_k_lake_recharge[per] - # zero-values to lak package lakes - data[per][model.isbc[0] == 1] = 0. - else: - if model.high_k_lake_recharge is not None: - # start with the last period with recharge data; update the high-k lake recharge - last_data_array[model.isbc[0] == 2] = model.high_k_lake_recharge[per] - last_data_array[model.isbc[0] == 1] = 0. - # assign to current per - data[per] = last_data_array + + # special case of LGR models + # with bottom connections to underlying parent cells + if model.version == 'mf6': + # (if model is an lgr inset model) + if model._is_lgr: + # regardless of what is specified for inset model bottom + # use top elevations of underlying parent model cells + nlay = model.cfg['dis']['dimensions']['nlay'] + lgr = model.parent.lgr[model.name] # Flopy Lgr inst. + n_refined = (np.array(lgr.ncppl) > 0).astype(int).sum() + # check if LGR model has connections to underlying parent cells + if (n_refined < model.parent.modelgrid.nlay): + # use the parent model bottoms + # mapped to the inset model grid by the Flopy Lgr util + #data[nlay-1] = lgr.botm[nlay-1] + # set the inset and parent model botms + # to the mean elevations of the inset cells + # connecting to each parent cell + ncpp = lgr.ncpp + nrowp = int(data[nlay-1].shape[0]/ncpp) + ncolp = int(data[nlay-1].shape[1]/ncpp) + lgr_inset_botm = data[nlay-1] + if any(np.isnan(lgr_inset_botm.flat)): + raise ValueError( + "Nan values in LGR inset model bottom; " + "can't use this to make the top of the parent model." + ) + # average each nccp x nccp block of inset model cells + topp = np.reshape(lgr_inset_botm, + (nrowp, ncpp, ncolp, ncpp)).mean(axis=(1, 3)) + n_parent_lgr_layers = np.sum(np.array(lgr.ncppl) > 0) + # remap averages back to the inset model shape + data[nlay-1] = np.repeat(np.repeat(topp, 5, axis=0), 5, axis=1) + # set the parent model botm in this area to be the same + lgr_area = model.parent.lgr[model.name].idomain == 0 + model.parent.dis.top[lgr_area[0]] = topp.ravel() + # set parent model layers in LGR area to zero-thickness + new_parent_botm = model.parent.dis.botm.array.copy() + for k in range(n_parent_lgr_layers): + new_parent_botm[k][lgr_area[0]] = topp.ravel() + model.parent.dis.botm = new_parent_botm + model.parent._update_top_botm_external_files() + + elif var in ['rech', 'recharge']: + if simulate_high_k_lakes: + for per in range(model.nper): + if per == 0 and per not in data: + raise KeyError("No recharge input specified for first stress period.") + if per in data: + # assign high-k lake recharge for stress period + # only assign if precip and open water evaporation data were read + # (otherwise keep original values in recharge array) + last_data_array = data[per].copy() + if model.high_k_lake_recharge is not None: + data[per][model.isbc[0] == 2] = model.high_k_lake_recharge[per] + # zero-values to lak package lakes + data[per][model.isbc[0] == 1] = 0. + else: + if model.high_k_lake_recharge is not None: + # start with the last period with recharge data; update the high-k lake recharge + last_data_array[model.isbc[0] == 2] = model.high_k_lake_recharge[per] + last_data_array[model.isbc[0] == 1] = 0. + # assign to current per + data[per] = last_data_array + if model.lgr: + for per, rech_array in data.items(): + rech_array[model._lgr_idomain2d == 0] = 0 elif var in ['hk', 'k'] and simulate_high_k_lakes: for i, arr in data.items(): data[i][model.isbc[i] == 2] = model.cfg['high_k_lakes']['high_k_value'] diff --git a/mfsetup/tests/data/pleasant_vertical_lgr_inset.yml b/mfsetup/tests/data/pleasant_vertical_lgr_inset.yml index ab163592..cb513b4e 100644 --- a/mfsetup/tests/data/pleasant_vertical_lgr_inset.yml +++ b/mfsetup/tests/data/pleasant_vertical_lgr_inset.yml @@ -1,10 +1,10 @@ simulation: - sim_name: 'pleasant_lgr' + sim_name: 'pleasant_vlgr' version: 'mf6' - sim_ws: 'pleasant_lgr/' + sim_ws: '../tmp/pleasant_vlgr/' model: - simulation: 'pleasant_lgr' + simulation: 'pleasant_vlgr' modelname: 'plsnt_lgr_inset' options: print_input: True @@ -32,7 +32,7 @@ model: setup_grid: source_data: features_shapefile: - filename: 'data/pleasant/source_data/shps/all_lakes.shp' + filename: '../../../examples/data/pleasant/source_data/shps/all_lakes.shp' id_column: 'HYDROID' include_ids: [600059060] # pleasant lake dxy: 40 # grid spacing, in CRS units (meters) @@ -44,16 +44,21 @@ dis: length_units: 'meters' dimensions: # nrow and ncol are based on the buffer distance and grid spacing - nlay: 5 + nlay: 9 source_data: top: - filename: 'data/pleasant/source_data/rasters/dem40m.tif' # DEM file; path relative to setup script + filename: '../../../examples/data/pleasant/source_data/rasters/dem40m.tif' # DEM file; path relative to setup script elevation_units: 'meters' botm: filenames: - 1: 'data/pleasant/source_data/rasters/botm0.tif' # preprocessed surface for parent model layer 0 bottom - 2: 'data/pleasant/source_data/rasters/botm1.tif' # preprocessed surface for parent model layer 1 bottom - 3: 'data/pleasant/source_data/rasters/botm2.tif' # preprocessed surface for parent model layer 2 bottom + # bottom of layer 0 set halfway between layer 1 bottom and model top (layer_refinement[0] == 2) + 1: '../../../examples/data/pleasant/source_data/rasters/botm0.tif' # preprocessed surface for parent model layer 0 bottom + # bottom of layer 2 set halfway between layers 1 and 3 (layer_refinement[1] == 2) + 3: '../../../examples/data/pleasant/source_data/rasters/botm1.tif' # preprocessed surface for parent model layer 1 bottom + # bottom of layer 4 set halfway between layers 3 and 5 (layer_refinement[2] == 2) + 5: '../../../examples/data/pleasant/source_data/rasters/botm2.tif' # preprocessed surface for parent model layer 2 bottom + # bottoms of layers 6 and 7 equally subdivide layers 5 and 8 (layer_refinement[3] == 3) + # bottom of layer 8 set from parent model layer 3 bottom # Recharge and Well packages are inherited from the parent model # (since this is an LGR model @@ -74,13 +79,13 @@ lak: # polygon shapefile of lake footprints littoral_zone_buffer_width: 40 lakes_shapefile: - filename: 'data/pleasant/source_data/shps/all_lakes.shp' + filename: '../../../examples/data/pleasant/source_data/shps/all_lakes.shp' id_column: 'HYDROID' include_ids: [600059060] # pleasant lake # setup lake precipitation and evaporation from PRISM data climate: filenames: - 600059060: 'data/pleasant/source_data/PRISM_ppt_tmean_stable_4km_189501_201901_43.9850_-89.5522.csv' + 600059060: '../../../examples/data/pleasant/source_data/PRISM_ppt_tmean_stable_4km_189501_201901_43.9850_-89.5522.csv' format: 'prism' period_stats: # for period 0, use average precip and evap for dates below @@ -90,11 +95,11 @@ lak: 1: 'mean' # average daily rate or value for each month # bathymetry file with lake depths to subtract off model top bathymetry_raster: - filename: 'data/pleasant/source_data/rasters/pleasant_bathymetry.tif' + filename: '../../../examples/data/pleasant/source_data/rasters/pleasant_bathymetry.tif' length_units: 'meters' # bathymetry file with stage/area/volume relationship stage_area_volume_file: - filename: 'data/pleasant/source_data/tables/area_stage_vol_Pleasant.csv' + filename: '../../../examples/data/pleasant/source_data/tables/area_stage_vol_Pleasant.csv' length_units: 'meters' id_column: 'hydroid' column_mappings: @@ -106,9 +111,9 @@ sfr: save_flows: True source_data: flowlines: - nhdplus_paths: ['data/pleasant/source_data/shps'] + nhdplus_paths: ['../../../examples/data/pleasant/source_data/shps'] dem: - filename: 'data/pleasant/source_data/rasters/dem40m.tif' + filename: '../../../examples/data/pleasant/source_data/rasters/dem40m.tif' elevation_units: 'meters' sfrmaker_options: set_streambed_top_elevations_from_dem: True @@ -118,11 +123,11 @@ obs: # any observations not within the LGR extent are culled source_data: filenames: [ - #'data/pleasant/source_data/tables/nwis_heads_info_file.csv', - 'data/pleasant/source_data/tables/lake_sites.csv', # some lakes are high-k; obs need to be in hydmod - 'data/pleasant/source_data/tables/wdnr_gw_sites.csv', - #'data/pleasant/source_data/tables/uwsp_heads.csv', - #'data/pleasant/source_data/tables/wgnhs_head_targets.csv' + #'../../../examples/data/pleasant/source_data/tables/nwis_heads_info_file.csv', + '../../../examples/data/pleasant/source_data/tables/lake_sites.csv', # some lakes are high-k; obs need to be in hydmod + '../../../examples/data/pleasant/source_data/tables/wdnr_gw_sites.csv', + #'../../../examples/data/pleasant/source_data/tables/uwsp_heads.csv', + #'../../../examples/data/pleasant/source_data/tables/wgnhs_head_targets.csv' ] column_mappings: obsname: ['obsprefix', 'common_name'] diff --git a/mfsetup/tests/data/pleasant_vertical_lgr_parent.yml b/mfsetup/tests/data/pleasant_vertical_lgr_parent.yml index e98514e8..bb234e4d 100644 --- a/mfsetup/tests/data/pleasant_vertical_lgr_parent.yml +++ b/mfsetup/tests/data/pleasant_vertical_lgr_parent.yml @@ -11,13 +11,13 @@ metadata: # input for MODFLOW 6 simulation simulation: - sim_name: 'pleasant_lgr' + sim_name: 'pleasant_vlgr' version: 'mf6' - sim_ws: 'pleasant_lgr/' + sim_ws: '../tmp/pleasant_vlgr/' # input for MODFLOW 6 model model: - simulation: 'pleasant_lgr' + simulation: 'pleasant_vlgr' modelname: 'plsnt_lgr_parent' options: print_input: True @@ -46,7 +46,7 @@ model: parent: # argments to flopy.modflow.Modflow.load for parent model namefile: 'pleasant.nam' - model_ws: 'data/pleasant/' + model_ws: '../../../examples/data/pleasant/' version: 'mfnwt' # information for modflow-setup # note: parent model is geo-located in namfile header @@ -56,11 +56,10 @@ parent: copy_stress_periods: 'all' inset_layer_mapping: # mapping between inset and parent model layers 0: 0 # inset: parent (inset layers 1 and 0 copied from parent layer 0) - 1: 0 - 2: 1 - 3: 2 + 1: 1 + 2: 2 + 3: 3 4: 3 - 5: 3 start_date_time: '2012-01-01' length_units: 'meters' time_units: 'days' @@ -79,14 +78,13 @@ setup_grid: # is set up in tandem and connected in the same MODFLOW 6 simulation lgr: pleasant_lgr_inset: - filename: 'pleasant_lgr_inset.yml' + filename: 'pleasant_vertical_lgr_inset.yml' # number of inset layers for each parent layer # this can be an integer (vertically uniform refinement) # a list (with an entry for each parent layer) - # or a dictionary (with entries for each parent layer that is refined; + # or a dictionary (with entries for each parent layer that is refined; # unlisted layers will default to zero) - # only a - layer_refinement: [1, 1, 1, 1, 1, 0] + layer_refinement: [2, 2, 2, 3, 0] # Structured Discretization Package dis: @@ -94,7 +92,7 @@ dis: length_units: 'meters' dimensions: # if nrow and ncol are not specified here, the entries above in setup_grid are used - nlay: 6 + nlay: 5 nrow: 25 ncol: 25 # griddata: @@ -103,14 +101,15 @@ dis: # idomain is developed from layer pinch-outs and boundary condition locations source_data: top: - filename: 'data/pleasant/source_data/rasters/dem40m.tif' # DEM file; path relative to setup script + filename: '../../../examples/data/pleasant/source_data/rasters/dem40m.tif' # DEM file; path relative to setup script elevation_units: 'meters' botm: filenames: - 1: 'data/pleasant/source_data/rasters/botm0.tif' # preprocessed surface for parent model layer 0 bottom - 2: 'data/pleasant/source_data/rasters/botm1.tif' # preprocessed surface for parent model layer 1 bottom - 3: 'data/pleasant/source_data/rasters/botm2.tif' # preprocessed surface for parent model layer 2 bottom - 5: 'data/pleasant/source_data/rasters/botm3.tif' # preprocessed surface for parent model layer 3 bottom + 0: '../../../examples/data/pleasant/source_data/rasters/botm0.tif' # preprocessed surface for parent model layer 0 bottom + 1: '../../../examples/data/pleasant/source_data/rasters/botm1.tif' # preprocessed surface for parent model layer 1 bottom + 2: '../../../examples/data/pleasant/source_data/rasters/botm2.tif' # preprocessed surface for parent model layer 2 bottom + # bottom of layer 3 set halfway between layers 2 and 5 + 4: '../../../examples/data/pleasant/source_data/rasters/botm3.tif' # preprocessed surface for parent model layer 3 bottom # Temporal Discretization Package tdis: @@ -141,7 +140,7 @@ ic: source_data: strt: from_parent: - binaryfile: 'data/pleasant/pleasant.hds' + binaryfile: '../../../examples/data/pleasant/pleasant.hds' stress_period: 0 # Node Property Flow Package @@ -180,16 +179,16 @@ sfr: source_data: flowlines: # path to NHDPlus version 2 dataset - nhdplus_paths: ['data/pleasant/source_data/shps'] + nhdplus_paths: ['../../../examples/data/pleasant/source_data/shps'] # if a DEM is included, streambed top elevations will be sampled # from the minimum DEM values within a 100 meter buffer around each stream reach dem: - filename: 'data/pleasant/source_data/rasters/dem40m.tif' + filename: '../../../examples/data/pleasant/source_data/rasters/dem40m.tif' elevation_units: 'meters' # SFR observations can be automatically setup from a CSV file # of x, y locations and fluxes observations: # see sfrmaker.observations.add_observations for arguments - filename: 'data/pleasant/source_data/tables/gages.csv' + filename: '../../../examples/data/pleasant/source_data/tables/gages.csv' obstype: 'downstream-flow' # modflow-6 observation type x_location_column: 'x' y_location_column: 'y' @@ -221,11 +220,11 @@ obs: # observations at each location can subsequently be processed # to represent a given screened interval # for example, using modflow-obs (https://github.com/aleaf/modflow-obs) - filenames: ['data/pleasant/source_data/tables/nwis_heads_info_file.csv', - 'data/pleasant/source_data/tables/lake_sites.csv', # some lakes are high-k; obs need to be in hydmod - 'data/pleasant/source_data/tables/wdnr_gw_sites.csv', - 'data/pleasant/source_data/tables/uwsp_heads.csv', - 'data/pleasant/source_data/tables/wgnhs_head_targets.csv' + filenames: ['../../../examples/data/pleasant/source_data/tables/nwis_heads_info_file.csv', + '../../../examples/data/pleasant/source_data/tables/lake_sites.csv', # some lakes are high-k; obs need to be in hydmod + '../../../examples/data/pleasant/source_data/tables/wdnr_gw_sites.csv', + '../../../examples/data/pleasant/source_data/tables/uwsp_heads.csv', + '../../../examples/data/pleasant/source_data/tables/wgnhs_head_targets.csv' ] column_mappings: obsname: ['obsprefix', 'obsnme', 'common_name'] @@ -235,4 +234,4 @@ obs: chd: perimeter_boundary: - parent_head_file: 'data/pleasant/pleasant.hds' # needed for the perimeter boundary setup + parent_head_file: '../../../examples/data/pleasant/pleasant.hds' # needed for the perimeter boundary setup diff --git a/mfsetup/tests/test_lgr.py b/mfsetup/tests/test_lgr.py index e8b0e64b..d1bc298b 100644 --- a/mfsetup/tests/test_lgr.py +++ b/mfsetup/tests/test_lgr.py @@ -1,11 +1,11 @@ import copy import glob import os -import shutil from pathlib import Path import flopy import numpy as np +import pandas as pd import pytest mf6 = flopy.mf6 @@ -15,7 +15,7 @@ from mfsetup import MF6model from mfsetup.discretization import make_lgr_idomain from mfsetup.fileio import dump, exe_exists, load, load_cfg, read_mf6_block -from mfsetup.interpolate import regrid3d +from mfsetup.grid import get_ij from mfsetup.mover import get_sfr_package_connections from mfsetup.testing import compare_inset_parent_values from mfsetup.utils import get_input_arguments @@ -25,6 +25,9 @@ def pleasant_lgr_test_cfg_path(project_root_path): return project_root_path + '/examples/pleasant_lgr_parent.yml' +@pytest.fixture(scope="session") +def pleasant_vertical_lgr_test_cfg_path(test_data_path): + return test_data_path / 'pleasant_vertical_lgr_parent.yml' @pytest.fixture(scope="function") def pleasant_lgr_cfg(pleasant_lgr_test_cfg_path): @@ -34,6 +37,13 @@ def pleasant_lgr_cfg(pleasant_lgr_test_cfg_path): cfg['gisdir'] = os.path.join(cfg['simulation']['sim_ws'], 'gis') return cfg +@pytest.fixture(scope="function") +def pleasant_vertical_lgr_cfg(pleasant_vertical_lgr_test_cfg_path): + cfg = load_cfg(pleasant_vertical_lgr_test_cfg_path, + default_file='/mf6_defaults.yml') + # add some stuff just for the tests + cfg['gisdir'] = os.path.join(cfg['simulation']['sim_ws'], 'gis') + return cfg @pytest.fixture(scope="function") def pleasant_simulation(pleasant_lgr_cfg): @@ -67,6 +77,11 @@ def pleasant_lgr_setup_from_yaml(pleasant_lgr_cfg): m.write_input() return m +@pytest.fixture(scope="function") +def pleasant_vertical_lgr_setup_from_yaml(pleasant_vertical_lgr_cfg): + m = MF6model.setup_from_cfg(pleasant_vertical_lgr_cfg) + m.write_input() + return m @pytest.fixture(scope="function") def pleasant_lgr_stand_alone_parent(pleasant_lgr_test_cfg_path, tmpdir): @@ -96,7 +111,8 @@ def pleasant_lgr_stand_alone_parent(pleasant_lgr_test_cfg_path, tmpdir): def test_make_lgr_idomain(get_pleasant_lgr_parent_with_grid): m = get_pleasant_lgr_parent_with_grid inset_model = m.inset['plsnt_lgr_inset'] - idomain = make_lgr_idomain(m.modelgrid, inset_model.modelgrid) + idomain = make_lgr_idomain(m.modelgrid, inset_model.modelgrid, + m.lgr['plsnt_lgr_inset'].ncppl) assert idomain.shape == m.modelgrid.shape l, b, r, t = inset_model.modelgrid.bounds isinset = (m.modelgrid.xcellcenters > l) & \ @@ -108,21 +124,34 @@ def test_make_lgr_idomain(get_pleasant_lgr_parent_with_grid): @pytest.mark.parametrize( - 'lgr_grid_spacing,parent_start_end_layers,inset_nlay', [ + 'lgr_grid_spacing,layer_refinement_input,inset_nlay', [ # parent layers 0 through 3 # are specified in pleasant_lgr_parent.yml - (40, 'use configuration', 4), - (40, (0, 2), 3), + (40, 'use configuration', 5), + (40, [1, 1, 1, 1, 0], 4), + (40, [2, 1, 1, 1, 0], 5), + # dictionary input option + (40, {0:1, 1:1, 2:1}, 3), + # integer layer refinement input option + (40, 1, 5), # special case to test setup - # with no parent start/end layers specified + # with no layer refinement specified (40, None, 5), pytest.param(35, None, 5, marks=pytest.mark.xfail( reason='inset spacing not a factor of parent spacing')), pytest.param(40, None, 4, marks=pytest.mark.xfail( reason='inset nlay inconsistent with parent layers')), + pytest.param(40, [1, 1, 1, 1], 5, marks=pytest.mark.xfail( + reason='List layer_refinement input must include a value for each layer')), + pytest.param(40, [0, 1, 1, 1, 1], 5, marks=pytest.mark.xfail( + reason='LGR child grid must start at model top')), + pytest.param(40, {1: 1}, 5, marks=pytest.mark.xfail( + reason='LGR child grid must start at model top')), + pytest.param(40, [1, 1, 0, 1, 1], 5, marks=pytest.mark.xfail( + reason='LGR child grid must be contiguous')), ] ) -def test_lgr_grid_setup(lgr_grid_spacing, parent_start_end_layers, +def test_lgr_grid_setup(lgr_grid_spacing, layer_refinement_input, inset_nlay, pleasant_lgr_cfg, pleasant_simulation, project_root_path): @@ -138,17 +167,14 @@ def test_lgr_grid_setup(lgr_grid_spacing, parent_start_end_layers, lgr_cfg = cfg['setup_grid']['lgr']['pleasant_lgr_inset'] lgr_cfg['cfg'] = inset_cfg del lgr_cfg['filename'] - if parent_start_end_layers is None: - del lgr_cfg['parent_start_layer'] - del lgr_cfg['parent_end_layer'] - k0, k1 = 0, None - elif parent_start_end_layers == 'use configuration': - k0 = lgr_cfg['parent_start_layer'] - k1 = lgr_cfg['parent_end_layer'] + if layer_refinement_input is None: + del lgr_cfg['layer_refinement'] + layer_refinement = [1] * cfg['dis']['dimensions']['nlay'] + elif layer_refinement_input == 'use configuration': + layer_refinement = lgr_cfg['layer_refinement'] else: - k0, k1 = parent_start_end_layers - lgr_cfg['parent_start_layer'] = k0 - lgr_cfg['parent_end_layer'] = k1 + lgr_cfg['layer_refinement'] = layer_refinement_input + layer_refinement = layer_refinement_input cfg['model']['simulation'] = pleasant_simulation kwargs = get_input_arguments(cfg['model'], mf6.ModflowGwf, exclude='packages') @@ -165,9 +191,14 @@ def test_lgr_grid_setup(lgr_grid_spacing, parent_start_end_layers, m.modelgrid.write_shapefile(outfolder / 'pleasant_lgr_parent_grid.shp') inset_model.modelgrid.write_shapefile(outfolder / 'pleasant_lgr_inset_grid.shp') + if np.isscalar(layer_refinement): + layer_refinement = np.array([1] * m.modelgrid.nlay) + elif isinstance(layer_refinement, dict): + layer_refinement = [layer_refinement.get(i, 0) for i in range(m.modelgrid.nlay)] + layer_refinement = np.array(layer_refinement) # verify that lgr area was removed from parent idomain lgr_idomain = make_lgr_idomain(m.modelgrid, inset_model.modelgrid, - k0, k1) + layer_refinement) idomain = m.idomain assert idomain[lgr_idomain == 0].sum() == 0 inset_cells_per_layer = inset_model.modelgrid.shape[1] *\ @@ -176,18 +207,28 @@ def test_lgr_grid_setup(lgr_grid_spacing, parent_start_end_layers, nparent_cells_within_lgr_per_layer = inset_cells_per_layer / refinement_factor**2 # for each layer where lgr is specified, # there should be at least enough inactive cells to cover the lrg grid area - if k1 is None: - k1 = m.modelgrid.nlay -1 layers_with_lgr = (lgr_idomain == 0).sum(axis=(1, 2)) >=\ nparent_cells_within_lgr_per_layer - assert all(layers_with_lgr[k0:k1]) + assert all(layers_with_lgr[layer_refinement > 0]) # layers outside of the specified lgr range should have # less inactive cells than the size of the lgr grid area # (specific to this test case with a large LGR extent relative to the total grid) - assert not any(layers_with_lgr[k1+1:]) - assert not any(layers_with_lgr[:k0]) - j=2 - # todo: add test that grids are aligned + assert not any(layers_with_lgr[layer_refinement == 0]) + + # make a plot for the docs + if layer_refinement_input == 'use configuration': + import matplotlib.pyplot as plt + fig, ax = plt.subplots(figsize=(10, 10)) + parent_mv = flopy.plot.PlotMapView(model=m, ax=ax, layer=0) + inset_mv = flopy.plot.PlotMapView(model=inset_model, ax=ax, layer=0) + l, r, b, t = m.modelgrid.extent + lcp = parent_mv.plot_grid(lw=0.5, ax=ax) + lci = inset_mv.plot_grid(lw=0.5) + ax.set_ylim(b, t) + ax.set_xlim(l, r) + ax.set_aspect(1) + plt.savefig(project_root_path + '/docs/source/_static/pleasant_lgr.png', + bbox_inches='tight') def test_setup_mover(pleasant_lgr_setup_from_yaml): @@ -201,6 +242,51 @@ def test_setup_mover(pleasant_lgr_setup_from_yaml): assert 'mover' in options +def test_lgr_parent_bcs_in_lgr_area(pleasant_vertical_lgr_setup_from_yaml): + """Test that boundary conditions specified in the parent model + inside of the LGR area don't get created + (these should be represented in the LGR inset/child model instead.""" + m = pleasant_vertical_lgr_setup_from_yaml + ghb_source_data = { + 'shapefile': { # pond near pleasant lake + 'filename': '../../../../examples/data/pleasant/source_data/shps/all_lakes.shp', + 'id_column': 'HYDROID', + 'include_ids': [600059161], + 'all_touched': True + }, + 'bhead': { + 'filename': '../../../../examples/data/pleasant/source_data/rasters/dem40m.tif', + 'elevation_units': 'meters', + 'stat': 'mean' + }, + 'cond': 9, # m2/d + } + m.setup_ghb(source_data=ghb_source_data) + # there should be no GHB Package, because this feature is inside of the LGR area + # (and should therefore be represented in the LGR inset/child model) + assert not hasattr(m, 'ghb') + wel_source_data = { + 'wells': { + # well with screen mostly in LGR area (should not be represented) + 'well1': {'per': 1, 'x': 555910.8, 'y': 389618.3, + 'screen_top': 300, 'screen_botm': 260, 'q': -2000}, + # well with screen mostly in parent model below LGR area (should be represented) + 'well2': {'per': 1, 'x': 555910.8, 'y': 389618.3, + 'screen_top': 250, 'screen_botm': 200, 'q': -2000} + } + } + added_wells = m.setup_wel(source_data=wel_source_data) + # well1 should not be in the parent model (open interval within LGR area) + assert 'well1' not in added_wells.stress_period_data.data[1]['boundname'] + assert 'well2' in added_wells.stress_period_data.data[1]['boundname'] + inset = m.inset['plsnt_lgr_inset'] + inset.setup_ghb(source_data=ghb_source_data) + assert hasattr(inset, 'ghb') + inset_added_wells = inset.setup_wel(source_data=wel_source_data) + assert 'well1' in inset_added_wells.stress_period_data.data[1]['boundname'] + assert 'well2' not in inset_added_wells.stress_period_data.data[1]['boundname'] + + def test_mover_get_sfr_package_connections(pleasant_lgr_setup_from_yaml): m = pleasant_lgr_setup_from_yaml parent_reach_data = m.sfrdata.reach_data @@ -225,6 +311,109 @@ def test_mover_get_sfr_package_connections(pleasant_lgr_setup_from_yaml): assert to_parent == {29: 8, 41: 1} +def test_lgr_bottom_elevations(pleasant_vertical_lgr_setup_from_yaml, mf6_exe): + """Test that boundary elevations specified in the LGR area + are the same as the top of the underlying active area in + the parent model (so there aren't gaps or overlap in the + numerical grid).""" + m = pleasant_vertical_lgr_setup_from_yaml + inset = m.inset['plsnt_lgr_inset'] + x = inset.modelgrid.xcellcenters + y = inset.modelgrid.ycellcenters + pi, pj = get_ij(m.modelgrid, x.ravel(), y.ravel(), + local=False) + pi = np.reshape(pi, x.shape) + pj = np.reshape(pj, x.shape) + parent_model_top = np.reshape(m.dis.top.array[pi, pj], x.shape) + + # LGR child bottom and parent top should be aligned + assert np.allclose(inset.dis.botm.array[-1], + parent_model_top) + + # check exchange data + exchange_data = pd.DataFrame(m.simulation.gwfgwf.exchangedata.array) + k1, i1, j1 = zip(*exchange_data['cellidm1']) + k2, i2, j2 = zip(*exchange_data['cellidm2']) + exchange_data = exchange_data.assign( + k1=k1, i1=i1, j1=j1, k2=k2, i2=i2, j2=j2) + vertical_connections = exchange_data.loc[exchange_data['ihc'] == 0] + + lgr = m.lgr[inset.name] # Flopy Lgr utility instance + active_cells = np.sum(inset.dis.idomain.array > 0, axis=0) > 0 + has_vertical_connection = np.zeros_like(active_cells) + has_vertical_connection[vertical_connections['i2'], vertical_connections['j2']] = 1 + # all active cells should have a vertical connection + assert all(has_vertical_connection[active_cells] == 1) + # each active horizontal location should only have 1 vertical connection + assert np.sum(has_vertical_connection[active_cells] == 1) == len(vertical_connections) + # check connection distances + all_layers = np.stack([m.dis.top.array] + [arr for arr in m.dis.botm.array]) + cell_thickness1 = -np.diff(all_layers, axis=0) + all_layers2 = np.stack([inset.dis.top.array] + [arr for arr in inset.dis.botm.array]) + cell_thickness2 = -np.diff(all_layers2, axis=0) + vertical_connections['layer_thick1'] = cell_thickness1[vertical_connections['k1'], + vertical_connections['i1'], + vertical_connections['j1']] + vertical_connections['layer_thick2'] = cell_thickness2[vertical_connections['k2'], + vertical_connections['i2'], + vertical_connections['j2']] + # cl1 should be 0.5 x parent cell thickness + assert np.allclose(vertical_connections['cl1'], vertical_connections['layer_thick1']/2) + # cl2 should be 0.5 x inset/child cell thickness + assert np.allclose(vertical_connections['cl2'], vertical_connections['layer_thick2']/2) + # hwva should equal inset cell bottom face area + inset_cell_areas = inset.dis.delc.array[vertical_connections['i2']] *\ + inset.dis.delr.array[vertical_connections['j2']] + assert np.allclose(vertical_connections['hwva'], inset_cell_areas) + + # test the horizontal connections too + horizontal_connections_ns = exchange_data.loc[(exchange_data['ihc'] > 0) & exchange_data['angldegx'].isin([0., 180.])] + horizontal_connections_ew = exchange_data.loc[(exchange_data['ihc'] > 0) & exchange_data['angldegx'].isin([90., 270.])] + # cl1 should be 0.5 x parent cell width + assert np.allclose(horizontal_connections_ns['cl1'], m.dis.delc.array[horizontal_connections_ns['i1']]/2) + assert np.allclose(horizontal_connections_ew['cl1'], m.dis.delr.array[horizontal_connections_ew['j1']]/2) + #cl2—is the distance between the center of cell 2 and the its shared face with cell 1. + assert np.allclose(horizontal_connections_ns['cl2'], inset.dis.delc.array[horizontal_connections_ns['i2']]/2) + assert np.allclose(horizontal_connections_ew['cl2'], inset.dis.delr.array[horizontal_connections_ew['j2']]/2) + #hwva is the parent cell width perpendicular to the connection + assert np.allclose(horizontal_connections_ns['hwva'], inset.dis.delr.array[horizontal_connections_ns['j2']]) + assert np.allclose(horizontal_connections_ew['hwva'], inset.dis.delc.array[horizontal_connections_ew['i2']]) + + # verify that the model runs + success = False + if exe_exists(mf6_exe): + success, buff = m.simulation.run_simulation() + if not success: + list_file = m.name_file.list.array + with open(list_file) as src: + list_output = src.read() + assert success, 'model run did not terminate successfully:\n{}'.format(list_output) + + # make a cross section plot for the documentation + import matplotlib.pyplot as plt + fig, ax = plt.subplots(figsize=(14, 5)) + # plot cross section through an area prone to mis-patch + # because of inconsistencies in input raster layer surfaces, + # and sub-200m variability in topography + #xs_line = [(553000, 390900), (558000, 390900)] + # line thru the lake + xs_line = [(553000, 390200), (558000, 390200)] + xs = flopy.plot.PlotCrossSection(model=m, + line={"line": xs_line}, ax=ax, + geographic_coords=True) + lc = xs.plot_grid(zorder=4) + xs2 = flopy.plot.PlotCrossSection(model=inset, + line={"line": xs_line}, ax=ax, + geographic_coords=True) + lc = xs2.plot_grid(zorder=4) + #xs2.plot_surface(inset.ic.strt.array, c='b') + #xs.plot_surface(m.ic.strt.array, c='b') + plt.savefig('../../../../docs/source/_static/pleasant_vlgr_xsection.png', + bbox_inches='tight') + plt.close() + + + def test_lgr_model_setup(pleasant_lgr_setup_from_yaml, tmpdir): m = pleasant_lgr_setup_from_yaml assert isinstance(m.inset, dict) diff --git a/mfsetup/tests/test_mf6_shellmound.py b/mfsetup/tests/test_mf6_shellmound.py index 768a0f17..6391d32d 100644 --- a/mfsetup/tests/test_mf6_shellmound.py +++ b/mfsetup/tests/test_mf6_shellmound.py @@ -349,11 +349,20 @@ def test_dis_setup(shellmound_model_with_grid): invalid_botms = np.ones_like(botm) invalid_botms[np.isnan(botm)] = 0 invalid_botms[thickness < 1.0001] = 0 + + import matplotlib.pyplot as plt + fig, axes = plt.subplots(3, 4) + for i, ax in enumerate(axes.flat): + if i < 12: + ax.imshow(m.idomain[i]) # these two arrays are not equal # because isolated cells haven't been removed from the second one # this verifies that _set_idomain is removing them - assert not np.array_equal(m.idomain[:, isactive].sum(axis=1), - invalid_botms[:, isactive].sum(axis=1)) + # note: 8/4/2024: after 012b6c6, this statement is True, + # presumably because the isolated cells are getting culled as part of setup_dis() + # TODO: might be better to test find_remove_isolated_cells directly + # assert not np.array_equal(m.idomain[:, isactive].sum(axis=1), + # invalid_botms[:, isactive].sum(axis=1)) invalid_botms = find_remove_isolated_cells(invalid_botms, minimum_cluster_size=20) active_cells = m.idomain[:, isactive].copy() active_cells[active_cells < 0] = 0 # need to do this because some idomain cells are -1 diff --git a/mfsetup/wells.py b/mfsetup/wells.py index 5be780bb..6226c51c 100644 --- a/mfsetup/wells.py +++ b/mfsetup/wells.py @@ -382,34 +382,6 @@ def assign_layers_from_screen_top_botm(data, model, T_r = T[::-1] data['k'] = len(T_r) - np.argmax(T_r, axis=0) - 1 - # get thicknesses for all layers - # (including portions of layers outside open interval) - all_layers = np.zeros((model.nlay + 1, model.nrow, model.ncol)) - all_layers[0] = model.dis.top.array - all_layers[1:] = model.dis.botm.array - all_layer_thicknesses = np.abs(np.diff(all_layers, axis=0)) - layer_thicknesses = -np.diff(all_layers[:, i, j], axis=0) - - # only include thicknesses for valid layers - # reset thicknesses to sat. thickness - if strt3D is not None: - sat_thickness = strt3D - model.dis.botm.array - # cells where the head is above the layer top - no_unsat = sat_thickness > all_layer_thicknesses - sat_thickness[no_unsat] = all_layer_thicknesses[no_unsat] - # cells where the head is below the cell bottom - sat_thickness[sat_thickness < 0] = 0 - layer_thicknesses = sat_thickness[:, i, j] - - # set inactive cells to 0 thickness for the purpose or relocating wells - layer_thicknesses[idomain[:, i, j] < 1] = 0 - data['idomain'] = idomain[data['k'], i, j] - data['laythick'] = layer_thicknesses[data['k'].values, - list(range(layer_thicknesses.shape[1]))] - # flag layers that are too thin or inactive - inactive = idomain[data.k, data.i, data.j] < 1 - invalid_open_interval = (data['laythick'] < minimum_layer_thickness) | inactive - outfile = model.cfg['wel']['output_files']['dropped_wells_file'].format(model.name) bad_wells = pd.DataFrame() # for LGR parent models, remove wells with >50% of their open interval within the LGR area @@ -438,6 +410,36 @@ def assign_layers_from_screen_top_botm(data, model, "should be checked)." ) data = data.loc[in_model].copy() + + + # get thicknesses for all layers + # (including portions of layers outside open interval) + k, i, j = data['k'].values, data['i'].values, data['j'].values + all_layers = np.zeros((model.nlay + 1, model.nrow, model.ncol)) + all_layers[0] = model.dis.top.array + all_layers[1:] = model.dis.botm.array + all_layer_thicknesses = np.abs(np.diff(all_layers, axis=0)) + layer_thicknesses = -np.diff(all_layers[:, i, j], axis=0) + + # only include thicknesses for valid layers + # reset thicknesses to sat. thickness + if strt3D is not None: + sat_thickness = strt3D - model.dis.botm.array + # cells where the head is above the layer top + no_unsat = sat_thickness > all_layer_thicknesses + sat_thickness[no_unsat] = all_layer_thicknesses[no_unsat] + # cells where the head is below the cell bottom + sat_thickness[sat_thickness < 0] = 0 + layer_thicknesses = sat_thickness[:, i, j] + + # set inactive cells to 0 thickness for the purpose or relocating wells + layer_thicknesses[idomain[:, i, j] < 1] = 0 + data['idomain'] = idomain[k, i, j] + data['laythick'] = layer_thicknesses[k, list(range(layer_thicknesses.shape[1]))] + # flag layers that are too thin or inactive + inactive = idomain[data.k, data.i, data.j] < 1 + invalid_open_interval = (data['laythick'] < minimum_layer_thickness) | inactive + if any(invalid_open_interval): # move wells that are still in a thin layer to the thickest active layer @@ -483,14 +485,14 @@ def assign_layers_from_screen_top_botm(data, model, minimum_layer_thickness, model.length_units)) if flux_col in data.columns: - pct_flux_below = 100 * bad_wells.loc[still_below_minimum, flux_col].sum()/data[flux_col].sum() + pct_flux_below = 100 * wells_in_too_thin_layers.loc[still_below_minimum, flux_col].sum()/data[flux_col].sum() msg += ', \nrepresenting {:.2f} %% of total flux,'.format(pct_flux_below) msg += '\nwere dropped. See {} for details.'.format(outfile) print(msg) # write shapefile and CSV output for wells that were dropped - bad_wells = pd.concat([bad_wells, data.loc[invalid_open_interval].copy()]) + bad_wells = pd.concat([bad_wells, wells_in_too_thin_layers]) bad_wells['routine'] = __name__ + '.assign_layers_from_screen_top_botm' cols = ['k', 'i', 'j', 'boundname', 'category', 'laythick', 'idomain', 'reason', 'routine', 'x', 'y']