From a07e8e54971466f861a80c25ae621fd8bfaa923b Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Thu, 25 Jul 2024 17:31:15 +0200 Subject: [PATCH 01/25] Start moving some logic to an abstract base class --- imod/mf6/wel.py | 279 +++++++++++++++++++++++++----------------------- 1 file changed, 143 insertions(+), 136 deletions(-) diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py index cc7cdb67f..6593ea16a 100644 --- a/imod/mf6/wel.py +++ b/imod/mf6/wel.py @@ -1,5 +1,6 @@ from __future__ import annotations +import abc import textwrap import warnings from datetime import datetime @@ -70,7 +71,148 @@ def mask_2D(package: Well, domain_2d: GridDataArray) -> Well: return cls._from_dataset(selection) -class Well(BoundaryCondition, IPointDataPackage): +class GridAgnosticWell(BoundaryCondition, IPointDataPackage, abc.ABC): + """ + Abstract base class for grid agnostic wells + """ + + @property + def x(self) -> npt.NDArray[np.float64]: + return self.dataset["x"].values + + @property + def y(self) -> npt.NDArray[np.float64]: + return self.dataset["y"].values + + @classmethod + def is_grid_agnostic_package(cls) -> bool: + return True + + def __create_cellid(self, wells_assigned: pd.DataFrame, active: xr.DataArray): + like = ones_like(active) + + # Groupby index and select first, to unset any duplicate records + # introduced by the multi-indexed "time" dimension. + df_for_cellid = wells_assigned.groupby("index").first() + d_for_cellid = df_for_cellid[["x", "y", "layer"]].to_dict("list") + + return self.__derive_cellid_from_points(like, **d_for_cellid) + + def __create_dataset_vars( + self, wells_assigned: pd.DataFrame, wells_df: pd.DataFrame, cellid: xr.DataArray + ) -> xr.Dataset: + """ + Create dataset with all variables (rate, concentration), with a similar shape as the cellids. + """ + data_vars = ["id", "rate"] + if "concentration" in wells_assigned.columns: + data_vars.append("concentration") + + ds_vars = wells_assigned[data_vars].to_xarray() + # "rate" variable in conversion from multi-indexed DataFrame to xarray + # DataArray results in duplicated values for "rate" along dimension + # "species". Select first species to reduce this again. + index_names = wells_df.index.names + if "species" in index_names: + ds_vars["rate"] = ds_vars["rate"].isel(species=0) + + # Carefully rename the dimension and set coordinates + d_rename = {"index": "ncellid"} + ds_vars = ds_vars.rename_dims(**d_rename).rename_vars(**d_rename) + ds_vars = ds_vars.assign_coords(**{"ncellid": cellid.coords["ncellid"].values}) + + return ds_vars + + @staticmethod + def __derive_cellid_from_points( + dst_grid: GridDataArray, + x: list, + y: list, + layer: list, + ) -> GridDataArray: + """ + Create DataArray with Modflow6 cell identifiers based on x, y coordinates + in a dataframe. For structured grid this DataArray contains 3 columns: + ``layer, row, column``. For unstructured grids, this contains 2 columns: + ``layer, cell2d``. + See also: https://water.usgs.gov/water-resources/software/MODFLOW-6/mf6io_6.4.0.pdf#page=35 + + Note + ---- + The "layer" coordinate should already be provided in the dataframe. + To determine the layer coordinate based on screen depts, look at + :func:`imod.prepare.wells.assign_wells`. + + Parameters + ---------- + dst_grid: {xr.DataArray, xu.UgridDataArray} + Destination grid to map the points to based on their x and y coordinates. + x: {list, np.array} + array-like with x-coordinates + y: {list, np.array} + array-like with y-coordinates + layer: {list, np.array} + array-like with layer-coordinates + + Returns + ------- + cellid : xr.DataArray + 2D DataArray with a ``ncellid`` rows and 3 to 2 columns, depending + on whether on a structured or unstructured grid.""" + + # Find indices belonging to x, y coordinates + indices_cell2d = points_indices(dst_grid, out_of_bounds="ignore", x=x, y=y) + # Convert cell2d indices from 0-based to 1-based. + indices_cell2d = {dim: index + 1 for dim, index in indices_cell2d.items()} + # Prepare layer indices, for later concatenation + + if isinstance(dst_grid, xu.UgridDataArray): + indices_layer = xr.DataArray( + layer, coords=indices_cell2d["mesh2d_nFaces"].coords + ) + face_dim = dst_grid.ugrid.grid.face_dimension + indices_cell2d_dims = [face_dim] + cell2d_coords = ["cell2d"] + else: + indices_layer = xr.DataArray(layer, coords=indices_cell2d["x"].coords) + indices_cell2d_dims = ["y", "x"] + cell2d_coords = ["row", "column"] + + # Prepare cellid array of the right shape. + cellid_ls = [indices_layer] + [ + indices_cell2d[dim] for dim in indices_cell2d_dims + ] + cellid = xr.concat(cellid_ls, dim="nmax_cellid") + # Rename generic dimension name "index" to ncellid. + cellid = cellid.rename(index="ncellid") + # Put dimensions in right order after concatenation. + cellid = cellid.transpose("ncellid", "nmax_cellid") + # Assign extra coordinate names. + coords = { + "nmax_cellid": ["layer"] + cell2d_coords, + "x": ("ncellid", x), + "y": ("ncellid", y), + } + cellid = cellid.assign_coords(**coords) + + return cellid + + def render(self, directory, pkgname, globaltimes, binary): + raise NotImplementedError( + f"{self.__class__.__name__} is a grid-agnostic package and does not have a render method. To render the package, first convert to a Modflow6 package by calling pkg.to_mf6_pkg()" + ) + + def write( + self, + pkgname: str, + globaltimes: Union[list[np.datetime64], np.ndarray], + write_context: WriteContext, + ): + raise NotImplementedError( + "To write a wel package first convert it to a MF6 well using to_mf6_pkg." + ) + +class Well(GridAgnosticWell): """ Agnostic WEL package, which accepts x, y and a top and bottom of the well screens. @@ -160,14 +302,6 @@ class Well(BoundaryCondition, IPointDataPackage): >>> imod.mf6.Well(x, y, screen_top, screen_bottom, rate_transient) """ - @property - def x(self) -> npt.NDArray[np.float64]: - return self.dataset["x"].values - - @property - def y(self) -> npt.NDArray[np.float64]: - return self.dataset["y"].values - _pkg_id = "wel" _auxiliary_data = {"concentration": "species"} @@ -242,10 +376,6 @@ def __init__( self.dataset = self.dataset.assign_coords(index=index_coord) self._validate_init_schemata(validate) - @classmethod - def is_grid_agnostic_package(cls) -> bool: - return True - def clip_box( self, time_min: Optional[cftime.datetime | np.datetime64 | str] = None, @@ -350,16 +480,6 @@ def _find_well_value_at_layer( return value - def write( - self, - pkgname: str, - globaltimes: Union[list[np.datetime64], np.ndarray], - write_context: WriteContext, - ): - raise NotImplementedError( - "To write a wel package first convert it to a MF6 well using to_mf6_pkg." - ) - def __create_wells_df(self) -> pd.DataFrame: wells_df = self.dataset.to_dataframe() wells_df = wells_df.rename( @@ -410,119 +530,6 @@ def __create_assigned_wells( return wells_assigned - def __create_dataset_vars( - self, wells_assigned: pd.DataFrame, wells_df: pd.DataFrame, cellid: xr.DataArray - ) -> xr.Dataset: - """ - Create dataset with all variables (rate, concentration), with a similar shape as the cellids. - """ - data_vars = ["id", "rate"] - if "concentration" in wells_assigned.columns: - data_vars.append("concentration") - - ds_vars = wells_assigned[data_vars].to_xarray() - # "rate" variable in conversion from multi-indexed DataFrame to xarray - # DataArray results in duplicated values for "rate" along dimension - # "species". Select first species to reduce this again. - index_names = wells_df.index.names - if "species" in index_names: - ds_vars["rate"] = ds_vars["rate"].isel(species=0) - - # Carefully rename the dimension and set coordinates - d_rename = {"index": "ncellid"} - ds_vars = ds_vars.rename_dims(**d_rename).rename_vars(**d_rename) - ds_vars = ds_vars.assign_coords(**{"ncellid": cellid.coords["ncellid"].values}) - - return ds_vars - - def __create_cellid(self, wells_assigned: pd.DataFrame, active: xr.DataArray): - like = ones_like(active) - - # Groupby index and select first, to unset any duplicate records - # introduced by the multi-indexed "time" dimension. - df_for_cellid = wells_assigned.groupby("index").first() - d_for_cellid = df_for_cellid[["x", "y", "layer"]].to_dict("list") - - return self.__derive_cellid_from_points(like, **d_for_cellid) - - @staticmethod - def __derive_cellid_from_points( - dst_grid: GridDataArray, - x: list, - y: list, - layer: list, - ) -> GridDataArray: - """ - Create DataArray with Modflow6 cell identifiers based on x, y coordinates - in a dataframe. For structured grid this DataArray contains 3 columns: - ``layer, row, column``. For unstructured grids, this contains 2 columns: - ``layer, cell2d``. - See also: https://water.usgs.gov/water-resources/software/MODFLOW-6/mf6io_6.4.0.pdf#page=35 - - Note - ---- - The "layer" coordinate should already be provided in the dataframe. - To determine the layer coordinate based on screen depts, look at - :func:`imod.prepare.wells.assign_wells`. - - Parameters - ---------- - dst_grid: {xr.DataArray, xu.UgridDataArray} - Destination grid to map the points to based on their x and y coordinates. - x: {list, np.array} - array-like with x-coordinates - y: {list, np.array} - array-like with y-coordinates - layer: {list, np.array} - array-like with layer-coordinates - - Returns - ------- - cellid : xr.DataArray - 2D DataArray with a ``ncellid`` rows and 3 to 2 columns, depending - on whether on a structured or unstructured grid.""" - - # Find indices belonging to x, y coordinates - indices_cell2d = points_indices(dst_grid, out_of_bounds="ignore", x=x, y=y) - # Convert cell2d indices from 0-based to 1-based. - indices_cell2d = {dim: index + 1 for dim, index in indices_cell2d.items()} - # Prepare layer indices, for later concatenation - - if isinstance(dst_grid, xu.UgridDataArray): - indices_layer = xr.DataArray( - layer, coords=indices_cell2d["mesh2d_nFaces"].coords - ) - face_dim = dst_grid.ugrid.grid.face_dimension - indices_cell2d_dims = [face_dim] - cell2d_coords = ["cell2d"] - else: - indices_layer = xr.DataArray(layer, coords=indices_cell2d["x"].coords) - indices_cell2d_dims = ["y", "x"] - cell2d_coords = ["row", "column"] - - # Prepare cellid array of the right shape. - cellid_ls = [indices_layer] + [ - indices_cell2d[dim] for dim in indices_cell2d_dims - ] - cellid = xr.concat(cellid_ls, dim="nmax_cellid") - # Rename generic dimension name "index" to ncellid. - cellid = cellid.rename(index="ncellid") - # Put dimensions in right order after concatenation. - cellid = cellid.transpose("ncellid", "nmax_cellid") - # Assign extra coordinate names. - coords = { - "nmax_cellid": ["layer"] + cell2d_coords, - "x": ("ncellid", x), - "y": ("ncellid", y), - } - cellid = cellid.assign_coords(**coords) - - return cellid - - def render(self, directory, pkgname, globaltimes, binary): - raise NotImplementedError( - f"{self.__class__.__name__} is a grid-agnostic package and does not have a render method. To render the package, first convert to a Modflow6 package by calling pkg.to_mf6_pkg()" - ) def to_mf6_pkg( self, From 4e0b8c682b359c4c4533ab72973c23ff45cd9933 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Fri, 26 Jul 2024 09:34:42 +0200 Subject: [PATCH 02/25] Make methods semi-private --- imod/mf6/wel.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py index 6593ea16a..e7a263377 100644 --- a/imod/mf6/wel.py +++ b/imod/mf6/wel.py @@ -88,7 +88,7 @@ def y(self) -> npt.NDArray[np.float64]: def is_grid_agnostic_package(cls) -> bool: return True - def __create_cellid(self, wells_assigned: pd.DataFrame, active: xr.DataArray): + def _create_cellid(self, wells_assigned: pd.DataFrame, active: xr.DataArray): like = ones_like(active) # Groupby index and select first, to unset any duplicate records @@ -96,9 +96,9 @@ def __create_cellid(self, wells_assigned: pd.DataFrame, active: xr.DataArray): df_for_cellid = wells_assigned.groupby("index").first() d_for_cellid = df_for_cellid[["x", "y", "layer"]].to_dict("list") - return self.__derive_cellid_from_points(like, **d_for_cellid) + return self._derive_cellid_from_points(like, **d_for_cellid) - def __create_dataset_vars( + def _create_dataset_vars( self, wells_assigned: pd.DataFrame, wells_df: pd.DataFrame, cellid: xr.DataArray ) -> xr.Dataset: """ @@ -124,7 +124,7 @@ def __create_dataset_vars( return ds_vars @staticmethod - def __derive_cellid_from_points( + def _derive_cellid_from_points( dst_grid: GridDataArray, x: list, y: list, @@ -603,9 +603,9 @@ def to_mf6_pkg( ) ds = xr.Dataset() - ds["cellid"] = self.__create_cellid(wells_assigned, active) + ds["cellid"] = self._create_cellid(wells_assigned, active) - ds_vars = self.__create_dataset_vars(wells_assigned, wells_df, ds["cellid"]) + ds_vars = self._create_dataset_vars(wells_assigned, wells_df, ds["cellid"]) ds = ds.assign(**ds_vars.data_vars) ds = remove_inactive(ds, active) From 75fcde6bc1b8668d9a098f57c6121419d07ed8be Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Fri, 26 Jul 2024 13:39:23 +0200 Subject: [PATCH 03/25] Move more methods to abstract well class --- imod/mf6/wel.py | 494 +++++++++++++++++++++++++++++++++++------------- 1 file changed, 364 insertions(+), 130 deletions(-) diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py index e7a263377..8a27917fe 100644 --- a/imod/mf6/wel.py +++ b/imod/mf6/wel.py @@ -61,7 +61,7 @@ def _assign_dims(arg: Any) -> Tuple | xr.DataArray: return "index", arg -def mask_2D(package: Well, domain_2d: GridDataArray) -> Well: +def mask_2D(package: GridAgnosticWell, domain_2d: GridDataArray) -> GridAgnosticWell: point_active = points_values(domain_2d, x=package.x, y=package.y) is_inside_exterior = point_active == 1 @@ -212,6 +212,138 @@ def write( "To write a wel package first convert it to a MF6 well using to_mf6_pkg." ) + def mask(self, domain: GridDataArray) -> GridAgnosticWell: + """ + Mask wells based on two-dimensional domain. For three-dimensional + masking: Wells falling in inactive cells are automatically removed in + the call to write to Modflow 6 package. You can verify this by calling + the ``to_mf6_pkg`` method. + """ + + # Drop layer coordinate if present, otherwise a layer coordinate is assigned + # which causes conflicts downstream when assigning wells and deriving + # cellids. + domain_2d = domain.isel(layer=0, drop=True, missing_dims="ignore").drop_vars( + "layer", errors="ignore" + ) + return mask_2D(self, domain_2d) + + def to_mf6_pkg( + self, + active: GridDataArray, + top: GridDataArray, + bottom: GridDataArray, + k: GridDataArray, + validate: bool = False, + is_partitioned: bool = False, + ) -> Mf6Wel: + """ + Write package to Modflow 6 package. + + Based on the model grid and top and bottoms, cellids are determined. + When well screens hit multiple layers, groundwater extractions are + distributed based on layer transmissivities. Wells located in inactive + cells are removed. + + Note + ---- + The well distribution based on transmissivities assumes confined + aquifers. If wells fall dry (and the rate distribution has to be + recomputed at runtime), it is better to use the Multi-Aquifer Well + package. + + Parameters + ---------- + active: {xarry.DataArray, xugrid.UgridDataArray} + Grid with active cells. + top: {xarry.DataArray, xugrid.UgridDataArray} + Grid with top of model layers. + bottom: {xarry.DataArray, xugrid.UgridDataArray} + Grid with bottom of model layers. + k: {xarry.DataArray, xugrid.UgridDataArray} + Grid with hydraulic conductivities. + validate: bool + Run validation before converting + is_partitioned: bool + Set to true if model has been partitioned + + Returns + ------- + Mf6Wel + Object with wells as list based input. + """ + if validate: + errors = self._validate(self._write_schemata) + if len(errors) > 0: + message = validation_pkg_error_message(errors) + raise ValidationError(message) + + wells_df = self._create_wells_df() + nwells_df = len(wells_df["id"].unique()) + wells_assigned = self._create_assigned_wells(wells_df, active, top, bottom, k) + + nwells_assigned = ( + 0 if wells_assigned.empty else len(wells_assigned["id"].unique()) + ) + + if nwells_df == 0: + raise ValueError("No wells were assigned in package. None were present.") + + if not is_partitioned and nwells_df != nwells_assigned: + raise ValueError( + "One or more well(s) are completely invalid due to minimum conductivity and thickness constraints." + ) + + ds = xr.Dataset() + ds["cellid"] = self._create_cellid(wells_assigned, active) + + ds_vars = self._create_dataset_vars(wells_assigned, wells_df, ds["cellid"]) + ds = ds.assign(**ds_vars.data_vars) + + ds = remove_inactive(ds, active) + ds["save_flows"] = self["save_flows"].values[()] + ds["print_flows"] = self["print_flows"].values[()] + ds["print_input"] = self["print_input"].values[()] + + filtered_wells = [ + id for id in wells_df["id"].unique() if id not in ds["id"].values + ] + if len(filtered_wells) > 0: + message = self.to_mf6_package_information(filtered_wells) + logger.log(loglevel=LogLevel.WARNING, message=message, additional_depth=2) + + ds = ds.drop_vars("id") + + return Mf6Wel(**ds.data_vars) + + def to_mf6_package_information(self, filtered_wells): + message = textwrap.dedent( + """Some wells were not placed in the MF6 well package. This + can be due to inactive cells or permeability/thickness constraints.\n""" + ) + if len(filtered_wells) < 10: + message += "The filtered wells are: \n" + else: + message += " The first 10 unplaced wells are: \n" + + for i in range(min(10, len(filtered_wells))): + message += f" id = {filtered_wells[i]} x = {self.dataset['x'][int(filtered_wells[i])].values[()]} y = {self.dataset['y'][int(filtered_wells[i])].values[()]} \n" + return message + + def _create_wells_df(self): + raise NotImplementedError("Method in abstract base class called") + + def _create_assigned_wells( + self, + wells_df: pd.DataFrame, + active: GridDataArray, + top: GridDataArray, + bottom: GridDataArray, + k: GridDataArray, + ): + raise NotImplementedError("Method in abstract base class called") + + class Well(GridAgnosticWell): """ Agnostic WEL package, which accepts x, y and a top and bottom of the well screens. @@ -419,7 +551,6 @@ def clip_box( y_max: optional, float top: optional, GridDataArray bottom: optional, GridDataArray - state_for_boundary: optional, GridDataArray Returns ------- @@ -480,7 +611,7 @@ def _find_well_value_at_layer( return value - def __create_wells_df(self) -> pd.DataFrame: + def _create_wells_df(self) -> pd.DataFrame: wells_df = self.dataset.to_dataframe() wells_df = wells_df.rename( columns={ @@ -496,15 +627,13 @@ def _validate(self, schemata: dict, **kwargs) -> dict[str, list[ValidationError] kwargs["screen_top"] = self.dataset["screen_top"] return Package._validate(self, schemata, **kwargs) - def __create_assigned_wells( + def _create_assigned_wells( self, wells_df: pd.DataFrame, active: GridDataArray, top: GridDataArray, bottom: GridDataArray, k: GridDataArray, - minimum_k: float, - minimum_thickness: float, ): # Ensure top, bottom & k # are broadcasted to 3d grid @@ -517,6 +646,9 @@ def __create_assigned_wells( index_names = wells_df.index.names + minimum_k = self.dataset["minimum_k"].item() + minimum_thickness = self.dataset["minimum_thickness"].item() + # Unset multi-index, because assign_wells cannot deal with # multi-indices which is returned by self.dataset.to_dataframe() in # case of a "time" and "species" coordinate. @@ -530,130 +662,6 @@ def __create_assigned_wells( return wells_assigned - - def to_mf6_pkg( - self, - active: GridDataArray, - top: GridDataArray, - bottom: GridDataArray, - k: GridDataArray, - validate: bool = False, - is_partitioned: bool = False, - ) -> Mf6Wel: - """ - Write package to Modflow 6 package. - - Based on the model grid and top and bottoms, cellids are determined. - When well screens hit multiple layers, groundwater extractions are - distributed based on layer transmissivities. Wells located in inactive - cells are removed. - - Note - ---- - The well distribution based on transmissivities assumes confined - aquifers. If wells fall dry (and the rate distribution has to be - recomputed at runtime), it is better to use the Multi-Aquifer Well - package. - - Parameters - ---------- - active: {xarry.DataArray, xugrid.UgridDataArray} - Grid with active cells. - top: {xarry.DataArray, xugrid.UgridDataArray} - Grid with top of model layers. - bottom: {xarry.DataArray, xugrid.UgridDataArray} - Grid with bottom of model layers. - k: {xarry.DataArray, xugrid.UgridDataArray} - Grid with hydraulic conductivities. - validate: bool - Run validation before converting - is_partitioned: bool - Set to true if model has been partitioned - - Returns - ------- - Mf6Wel - Object with wells as list based input. - """ - if validate: - errors = self._validate(self._write_schemata) - if len(errors) > 0: - message = validation_pkg_error_message(errors) - raise ValidationError(message) - - minimum_k = self.dataset["minimum_k"].item() - minimum_thickness = self.dataset["minimum_thickness"].item() - - wells_df = self.__create_wells_df() - nwells_df = len(wells_df["id"].unique()) - wells_assigned = self.__create_assigned_wells( - wells_df, active, top, bottom, k, minimum_k, minimum_thickness - ) - - nwells_assigned = ( - 0 if wells_assigned.empty else len(wells_assigned["id"].unique()) - ) - - if nwells_df == 0: - raise ValueError("No wells were assigned in package. None were present.") - - if not is_partitioned and nwells_df != nwells_assigned: - raise ValueError( - "One or more well(s) are completely invalid due to minimum conductivity and thickness constraints." - ) - - ds = xr.Dataset() - ds["cellid"] = self._create_cellid(wells_assigned, active) - - ds_vars = self._create_dataset_vars(wells_assigned, wells_df, ds["cellid"]) - ds = ds.assign(**ds_vars.data_vars) - - ds = remove_inactive(ds, active) - ds["save_flows"] = self["save_flows"].values[()] - ds["print_flows"] = self["print_flows"].values[()] - ds["print_input"] = self["print_input"].values[()] - - filtered_wells = [ - id for id in wells_df["id"].unique() if id not in ds["id"].values - ] - if len(filtered_wells) > 0: - message = self.to_mf6_package_information(filtered_wells) - logger.log(loglevel=LogLevel.WARNING, message=message, additional_depth=2) - - ds = ds.drop_vars("id") - - return Mf6Wel(**ds.data_vars) - - def to_mf6_package_information(self, filtered_wells): - message = textwrap.dedent( - """Some wells were not placed in the MF6 well package. This - can be due to inactive cells or permeability/thickness constraints.\n""" - ) - if len(filtered_wells) < 10: - message += "The filtered wells are: \n" - else: - message += " The first 10 unplaced wells are: \n" - - for i in range(min(10, len(filtered_wells))): - message += f" id = {filtered_wells[i]} x = {self.dataset['x'][int(filtered_wells[i])].values[()]} y = {self.dataset['y'][int(filtered_wells[i])].values[()]} \n" - return message - - def mask(self, domain: GridDataArray) -> Well: - """ - Mask wells based on two-dimensional domain. For three-dimensional - masking: Wells falling in inactive cells are automatically removed in - the call to write to Modflow 6 package. You can verify this by calling - the ``to_mf6_pkg`` method. - """ - - # Drop layer coordinate if present, otherwise a layer coordinate is assigned - # which causes conflicts downstream when assigning wells and deriving - # cellids. - domain_2d = domain.isel(layer=0, drop=True, missing_dims="ignore").drop_vars( - "layer", errors="ignore" - ) - return mask_2D(self, domain_2d) - @classmethod def from_imod5_data( cls, @@ -731,6 +739,232 @@ def from_imod5_data( ) +class LayeredWell(GridAgnosticWell): + """ + Agnostic WEL package, which accepts x, y and layers. + + This package can be written to any provided model grid. + Any number of WEL Packages can be specified for a single groundwater flow model. + https://water.usgs.gov/water-resources/software/MODFLOW-6/mf6io_6.0.4.pdf#page=63 + + Parameters + ---------- + + y: float or list of floats or np.array of floats + is the y location of the well. + x: float or list of floats or np.array of floats + is the x location of the well. + layer: int or list of ints or np.array of ints + is the layer of the well. + rate: float, list of floats or xr.DataArray + is the volumetric well rate. A positive value indicates well + (injection) and a negative value indicates discharge (extraction) (q). + If provided as DataArray, an ``"index"`` dimension is required and an + optional ``"time"`` dimension and coordinate specify transient input. + In the latter case, it is important that dimensions are in the order: + ``("time", "index")`` + concentration: array of floats (xr.DataArray, optional) + if this flow package is used in simulations also involving transport, then this array is used + as the concentration for inflow over this boundary. + concentration_boundary_type: ({"AUX", "AUXMIXED"}, optional) + if this flow package is used in simulations also involving transport, then this keyword specifies + how outflow over this boundary is computed. + id: list of Any, optional + assign an identifier code to each well. if not provided, one will be generated + Must be convertible to string, and unique entries. + minimum_k: float, optional + on creating point wells, no point wells will be placed in cells with a lower horizontal conductivity than this + minimum_thickness: float, optional + on creating point wells, no point wells will be placed in cells with a lower thickness than this + print_input: ({True, False}, optional) + keyword to indicate that the list of well information will be written to + the listing file immediately after it is read. + Default is False. + print_flows: ({True, False}, optional) + Indicates that the list of well flow rates will be printed to the + listing file for every stress period time step in which "BUDGET PRINT" + is specified in Output Control. If there is no Output Control option + and PRINT FLOWS is specified, then flow rates are printed for the last + time step of each stress period. + Default is False. + save_flows: ({True, False}, optional) + Indicates that well flow terms will be written to the file specified + with "BUDGET FILEOUT" in Output Control. + Default is False. + observations: [Not yet supported.] + Default is None. + validate: {True, False} + Flag to indicate whether the package should be validated upon + initialization. This raises a ValidationError if package input is + provided in the wrong manner. Defaults to True. + repeat_stress: Optional[xr.DataArray] of datetimes + Used to repeat data for e.g. repeating stress periods such as + seasonality without duplicating the values. The DataArray should have + dimensions ``("repeat", "repeat_items")``. The ``repeat_items`` + dimension should have size 2: the first value is the "key", the second + value is the "value". For the "key" datetime, the data of the "value" + datetime will be used. Can also be set with a dictionary using the + ``set_repeat_stress`` method. + + Examples + --------- + + >>> layer = [1, 2] + >>> y = [83.0, 77.0] + >>> x = [81.0, 82.0] + >>> rate = [1.0, 1.0] + + >>> imod.mf6.Well(x, y, layer, rate) + + For a transient well: + + >>> weltimes = pd.date_range("2000-01-01", "2000-01-03") + + >>> rate_factor_time = xr.DataArray([0.5, 1.0], coords={"time": weltimes}, dims=("time",)) + >>> rate_transient = rate_factor_time * xr.DataArray(rate, dims=("index",)) + + >>> imod.mf6.Well(x, y, layer, rate_transient) + """ + + _pkg_id = "wel" + + _auxiliary_data = {"concentration": "species"} + _init_schemata = { + "layer": [DTypeSchema(np.integer)], + "y": [DTypeSchema(np.floating)], + "x": [DTypeSchema(np.floating)], + "rate": [DTypeSchema(np.floating)], + "concentration": [DTypeSchema(np.floating)], + } + _write_schemata = { + "layer": [AnyNoDataSchema(), EmptyIndexesSchema()], + "y": [AnyNoDataSchema(), EmptyIndexesSchema()], + "x": [AnyNoDataSchema(), EmptyIndexesSchema()], + "rate": [AnyNoDataSchema(), EmptyIndexesSchema()], + "concentration": [AnyNoDataSchema(), EmptyIndexesSchema()], + } + + @init_log_decorator() + def __init__( + self, + x: np.ndarray | list[float], + y: np.ndarray | list[float], + layer: np.ndarray | list[int], + rate: list[float] | xr.DataArray, + concentration: Optional[list[float] | xr.DataArray] = None, + concentration_boundary_type="aux", + id: Optional[list[Any]] = None, + minimum_k: float = 0.1, + minimum_thickness: float = 1.0, + print_input: bool = False, + print_flows: bool = False, + save_flows: bool = False, + observations=None, + validate: bool = True, + repeat_stress: Optional[xr.DataArray] = None, + ): + if id is None: + id = [str(i) for i in range(len(x))] + else: + set_id = set(id) + if len(id) != len(set_id): + raise ValueError("id's must be unique") + id = [str(i) for i in id] + dict_dataset = { + "layer": _assign_dims(layer), + "y": _assign_dims(y), + "x": _assign_dims(x), + "rate": _assign_dims(rate), + "id": _assign_dims(id), + "minimum_k": minimum_k, + "minimum_thickness": minimum_thickness, + "print_input": print_input, + "print_flows": print_flows, + "save_flows": save_flows, + "observations": observations, + "repeat_stress": repeat_stress, + "concentration": concentration, + "concentration_boundary_type": concentration_boundary_type, + } + super().__init__(dict_dataset) + # Set index as coordinate + index_coord = np.arange(self.dataset.dims["index"]) + self.dataset = self.dataset.assign_coords(index=index_coord) + self._validate_init_schemata(validate) + + def clip_box( + self, + time_min: Optional[cftime.datetime | np.datetime64 | str] = None, + time_max: Optional[cftime.datetime | np.datetime64 | str] = None, + layer_min: Optional[int] = None, + layer_max: Optional[int] = None, + x_min: Optional[float] = None, + x_max: Optional[float] = None, + y_min: Optional[float] = None, + y_max: Optional[float] = None, + top: Optional[GridDataArray] = None, + bottom: Optional[GridDataArray] = None, + ) -> Package: + """ + Clip a package by a bounding box (time, layer, y, x). + + Slicing intervals may be half-bounded, by providing None: + + * To select 500.0 <= x <= 1000.0: + ``clip_box(x_min=500.0, x_max=1000.0)``. + * To select x <= 1000.0: ``clip_box(x_min=None, x_max=1000.0)`` + or ``clip_box(x_max=1000.0)``. + * To select x >= 500.0: ``clip_box(x_min = 500.0, x_max=None.0)`` + or ``clip_box(x_min=1000.0)``. + + Parameters + ---------- + time_min: optional + time_max: optional + layer_min: optional, int + layer_max: optional, int + x_min: optional, float + x_max: optional, float + y_min: optional, float + y_max: optional, float + top: optional, GridDataArray + bottom: optional, GridDataArray + + Returns + ------- + sliced : Package + """ + # The super method will select in the time dimension without issues. + new = super().clip_box(time_min=time_min, time_max=time_max) + + ds = new.dataset + + # Initiate array of True with right shape to deal with case no spatial + # selection needs to be done. + in_bounds = np.full(ds.dims["index"], True) + # Select all variables along "index" dimension + in_bounds &= values_within_range(ds["x"], x_min, x_max) + in_bounds &= values_within_range(ds["y"], y_min, y_max) + in_bounds &= values_within_range(ds["layer"], layer_min, layer_max) + # Replace dataset with reduced dataset based on booleans + new.dataset = ds.loc[{"index": in_bounds}] + + return new + + def _create_wells_df(self) -> pd.DataFrame: + return self.dataset.to_dataframe() + + def _create_assigned_wells( + self, + wells_df: pd.DataFrame, + active: GridDataArray, + top: GridDataArray, + bottom: GridDataArray, + k: GridDataArray, + ): + raise NotImplementedError("Method in abstract base class called") + + class WellDisStructured(DisStructuredBoundaryCondition): """ WEL package for structured discretization (DIS) models . From 58512fca58031159f12d5e8204c3fa907f9b235a Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Fri, 26 Jul 2024 13:41:16 +0200 Subject: [PATCH 04/25] Format code a bit to avoid long line format string --- imod/mf6/wel.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py index 8a27917fe..1619a695e 100644 --- a/imod/mf6/wel.py +++ b/imod/mf6/wel.py @@ -327,7 +327,10 @@ def to_mf6_package_information(self, filtered_wells): message += " The first 10 unplaced wells are: \n" for i in range(min(10, len(filtered_wells))): - message += f" id = {filtered_wells[i]} x = {self.dataset['x'][int(filtered_wells[i])].values[()]} y = {self.dataset['y'][int(filtered_wells[i])].values[()]} \n" + ids = filtered_wells[i] + x = self.dataset['x'][int(filtered_wells[i])].values[()] + y = self.dataset['y'][int(filtered_wells[i])].values[()] + message += f" id = {ids} x = {x} y = {y} \n" return message def _create_wells_df(self): From f037f58b603f20f6f65b518702b4dbdf8a849505 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Fri, 26 Jul 2024 14:21:29 +0200 Subject: [PATCH 05/25] First version of finished LayeredWell --- imod/mf6/wel.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py index 1619a695e..7c0fd9e82 100644 --- a/imod/mf6/wel.py +++ b/imod/mf6/wel.py @@ -328,8 +328,8 @@ def to_mf6_package_information(self, filtered_wells): for i in range(min(10, len(filtered_wells))): ids = filtered_wells[i] - x = self.dataset['x'][int(filtered_wells[i])].values[()] - y = self.dataset['y'][int(filtered_wells[i])].values[()] + x = self.dataset["x"][int(filtered_wells[i])].values[()] + y = self.dataset["y"][int(filtered_wells[i])].values[()] message += f" id = {ids} x = {x} y = {y} \n" return message @@ -644,7 +644,6 @@ def _create_assigned_wells( bottom = like * bottom top_2d = (like * top).sel(layer=1) top_3d = bottom.shift(layer=1).fillna(top_2d) - k = like * k index_names = wells_df.index.names @@ -965,7 +964,7 @@ def _create_assigned_wells( bottom: GridDataArray, k: GridDataArray, ): - raise NotImplementedError("Method in abstract base class called") + return wells_df class WellDisStructured(DisStructuredBoundaryCondition): From a188d4882c2ad87739757fb35b80a83384841ce6 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Fri, 26 Jul 2024 16:16:32 +0200 Subject: [PATCH 06/25] Refactor code into separate cases and add layered well cases and tests. --- imod/mf6/__init__.py | 2 +- imod/tests/test_mf6/test_mf6_wel.py | 314 +++++++++++++++++----------- 2 files changed, 189 insertions(+), 127 deletions(-) diff --git a/imod/mf6/__init__.py b/imod/mf6/__init__.py index d6225b065..652602432 100644 --- a/imod/mf6/__init__.py +++ b/imod/mf6/__init__.py @@ -49,5 +49,5 @@ from imod.mf6.timedis import TimeDiscretization from imod.mf6.utilities.regrid import RegridderType, RegridderWeightsCache from imod.mf6.uzf import UnsaturatedZoneFlow -from imod.mf6.wel import Well, WellDisStructured, WellDisVertices +from imod.mf6.wel import LayeredWell, Well, WellDisStructured, WellDisVertices from imod.mf6.write_context import WriteContext diff --git a/imod/tests/test_mf6/test_mf6_wel.py b/imod/tests/test_mf6/test_mf6_wel.py index 50deea31b..85f1057cc 100644 --- a/imod/tests/test_mf6/test_mf6_wel.py +++ b/imod/tests/test_mf6/test_mf6_wel.py @@ -34,39 +34,154 @@ ] -def test_to_mf6_pkg__high_lvl_stationary(basic_dis, well_high_lvl_test_data_stationary): +class GridAgnosticWellCases: + def case_well_stationary(self, well_high_lvl_test_data_stationary): + obj = imod.mf6.Well(*well_high_lvl_test_data_stationary) + dims_expected = { + "ncellid": 8, + "nmax_cellid": 3, + "species": 2, + } + cellid_expected = np.array( + [ + [1, 1, 9], + [1, 2, 9], + [1, 1, 8], + [1, 2, 8], + [2, 3, 7], + [2, 4, 7], + [2, 3, 6], + [2, 4, 6], + ], + dtype=np.int64, + ) + rate_expected = np.array(np.ones((8,), dtype=np.float32)) + return obj, dims_expected, cellid_expected, rate_expected + + def case_well_stationary_multilevel(self, well_high_lvl_test_data_stationary): + x, y, screen_top, _, rate_wel, concentration = well_high_lvl_test_data_stationary + screen_bottom = [-20.0] * 8 + obj = imod.mf6.Well(x, y, screen_top, screen_bottom, rate_wel, concentration) + dims_expected = { + "ncellid": 12, + "nmax_cellid": 3, + "species": 2, + } + cellid_expected = np.array( + [ + [1, 1, 9], + [1, 2, 9], + [1, 1, 8], + [1, 2, 8], + [2, 1, 9], + [2, 2, 9], + [2, 1, 8], + [2, 2, 8], + [2, 3, 7], + [2, 4, 7], + [2, 3, 6], + [2, 4, 6], + ], + dtype=np.int64, + ) + rate_expected = np.array([0.25] * 4 + [0.75] * 4 + [1.0] * 4) + return obj, dims_expected, cellid_expected, rate_expected + + def case_well_transient(self, well_high_lvl_test_data_transient): + obj = imod.mf6.Well(*well_high_lvl_test_data_transient) + dims_expected = { + "ncellid": 8, + "time": 5, + "nmax_cellid": 3, + "species": 2, + } + cellid_expected = np.array( + [ + [1, 1, 9], + [1, 2, 9], + [1, 1, 8], + [1, 2, 8], + [2, 3, 7], + [2, 4, 7], + [2, 3, 6], + [2, 4, 6], + ], + dtype=np.int64, + ) + rate_expected = np.outer(np.ones((8,), dtype=np.float32), np.arange(5) + 1) + return obj, dims_expected, cellid_expected, rate_expected + + def case_layered_well_stationary(self, well_high_lvl_test_data_stationary): + x, y, _, _, rate_wel, concentration = well_high_lvl_test_data_stationary + layer = [1, 1, 1, 1, 2, 2, 2, 2] + obj = imod.mf6.LayeredWell(x, y, layer, rate_wel, concentration) + dims_expected = { + "ncellid": 8, + "nmax_cellid": 3, + "species": 2, + } + cellid_expected = np.array( + [ + [1, 1, 9], + [1, 2, 9], + [1, 1, 8], + [1, 2, 8], + [2, 3, 7], + [2, 4, 7], + [2, 3, 6], + [2, 4, 6], + ], + dtype=np.int64, + ) + rate_expected = np.array(np.ones((8,), dtype=np.float32)) + return obj, dims_expected, cellid_expected, rate_expected + + def case_layered_well_transient(self, well_high_lvl_test_data_transient): + x, y, _, _, rate_wel, concentration = well_high_lvl_test_data_transient + layer = [1, 1, 1, 1, 2, 2, 2, 2] + obj = imod.mf6.LayeredWell(x, y, layer, rate_wel, concentration) + dims_expected = { + "ncellid": 8, + "time": 5, + "nmax_cellid": 3, + "species": 2, + } + cellid_expected = np.array( + [ + [1, 1, 9], + [1, 2, 9], + [1, 1, 8], + [1, 2, 8], + [2, 3, 7], + [2, 4, 7], + [2, 3, 6], + [2, 4, 6], + ], + dtype=np.int64, + ) + rate_expected = np.outer(np.ones((8,), dtype=np.float32), np.arange(5) + 1) + return obj, dims_expected, cellid_expected, rate_expected + + +@parametrize_with_cases( + ["wel", "dims_expected", "cellid_expected", "rate_expected"], cases=GridAgnosticWellCases +) +def test_to_mf6_pkg( + basic_dis, wel, dims_expected, cellid_expected, rate_expected +): # Arrange idomain, top, bottom = basic_dis - wel = imod.mf6.Well(*well_high_lvl_test_data_stationary) active = idomain == 1 k = xr.ones_like(idomain) nmax_cellid_expected = np.array(["layer", "row", "column"]) - cellid_expected = np.array( - [ - [1, 1, 9], - [1, 2, 9], - [1, 1, 8], - [1, 2, 8], - [2, 3, 7], - [2, 4, 7], - [2, 3, 6], - [2, 4, 6], - ], - dtype=np.int64, - ) - rate_expected = np.array(np.ones((8,), dtype=np.float32)) # Act mf6_wel = wel.to_mf6_pkg(active, top, bottom, k) mf6_ds = mf6_wel.dataset # Assert - assert dict(mf6_ds.dims) == { - "ncellid": 8, - "nmax_cellid": 3, - "species": 2, - } + assert dict(mf6_ds.dims) == dims_expected np.testing.assert_equal(mf6_ds.coords["nmax_cellid"].values, nmax_cellid_expected) np.testing.assert_equal(mf6_ds["cellid"].values, cellid_expected) np.testing.assert_equal(mf6_ds["rate"].values, rate_expected) @@ -109,94 +224,6 @@ def test_to_mf6_pkg__validate_filter_top(well_high_lvl_test_data_stationary): ) -def test_to_mf6_pkg__high_lvl_multilevel(basic_dis, well_high_lvl_test_data_stationary): - """ - Test with stationary wells where the first 4 well screens extend over 2 layers. - Rates are distributed based on the fraction of the screen length in each layer. - In this case: The first layer should get 0.25, the second 0.75. - """ - # Arrange - idomain, top, bottom = basic_dis - x, y, screen_top, _, rate_wel, concentration = well_high_lvl_test_data_stationary - screen_bottom = [-20.0] * 8 - wel = imod.mf6.Well(x, y, screen_top, screen_bottom, rate_wel, concentration) - active = idomain == 1 - k = xr.ones_like(idomain) - - nmax_cellid_expected = np.array(["layer", "row", "column"]) - cellid_expected = np.array( - [ - [1, 1, 9], - [1, 2, 9], - [1, 1, 8], - [1, 2, 8], - [2, 1, 9], - [2, 2, 9], - [2, 1, 8], - [2, 2, 8], - [2, 3, 7], - [2, 4, 7], - [2, 3, 6], - [2, 4, 6], - ], - dtype=np.int64, - ) - rate_expected = np.array([0.25] * 4 + [0.75] * 4 + [1.0] * 4) - - # Act - mf6_wel = wel.to_mf6_pkg(active, top, bottom, k) - mf6_ds = mf6_wel.dataset - - # Assert - assert dict(mf6_ds.dims) == { - "ncellid": 12, - "nmax_cellid": 3, - "species": 2, - } - np.testing.assert_equal(mf6_ds.coords["nmax_cellid"].values, nmax_cellid_expected) - np.testing.assert_equal(mf6_ds["cellid"].values, cellid_expected) - np.testing.assert_equal(mf6_ds["rate"].values, rate_expected) - - -def test_to_mf6_pkg__high_lvl_transient(basic_dis, well_high_lvl_test_data_transient): - # Arrange - idomain, top, bottom = basic_dis - wel = imod.mf6.Well(*well_high_lvl_test_data_transient) - active = idomain == 1 - k = xr.ones_like(idomain) - - nmax_cellid_expected = np.array(["layer", "row", "column"]) - cellid_expected = np.array( - [ - [1, 1, 9], - [1, 2, 9], - [1, 1, 8], - [1, 2, 8], - [2, 3, 7], - [2, 4, 7], - [2, 3, 6], - [2, 4, 6], - ], - dtype=np.int64, - ) - rate_expected = np.outer(np.ones((8,), dtype=np.float32), np.arange(5) + 1) - - # Act - mf6_wel = wel.to_mf6_pkg(active, top, bottom, k) - mf6_ds = mf6_wel.dataset - - # Assert - assert dict(mf6_wel.dataset.dims) == { - "ncellid": 8, - "time": 5, - "nmax_cellid": 3, - "species": 2, - } - np.testing.assert_equal(mf6_ds.coords["nmax_cellid"].values, nmax_cellid_expected) - np.testing.assert_equal(mf6_ds["cellid"].values, cellid_expected) - np.testing.assert_equal(mf6_ds["rate"].values, rate_expected) - - def test_to_mf6_pkg__logging_with_message( tmp_path, basic_dis, well_high_lvl_test_data_transient ): @@ -350,7 +377,7 @@ def case_clip_xy(parameterizable_basic_dis): } expected_dims = {"index": 3, "species": 2} - return clip_arguments, expected_dims, does_not_raise() + return clip_arguments, expected_dims, does_not_raise(), does_not_raise() @staticmethod def case_clip_layer_max(parameterizable_basic_dis): @@ -358,7 +385,7 @@ def case_clip_layer_max(parameterizable_basic_dis): clip_arguments = {"layer_max": 2, "bottom": bottom, "top": top} expected_dims = {"index": 4, "species": 2} - return clip_arguments, expected_dims, does_not_raise() + return clip_arguments, expected_dims, does_not_raise(), does_not_raise() @staticmethod def case_clip_layer_min(parameterizable_basic_dis): @@ -366,7 +393,7 @@ def case_clip_layer_min(parameterizable_basic_dis): clip_arguments = {"layer_min": 5, "bottom": bottom, "top": top} expected_dims = {"index": 4, "species": 2} - return clip_arguments, expected_dims, does_not_raise() + return clip_arguments, expected_dims, does_not_raise(), does_not_raise() @staticmethod def case_clip_layer_min_layer_max(parameterizable_basic_dis): @@ -374,7 +401,7 @@ def case_clip_layer_min_layer_max(parameterizable_basic_dis): clip_arguments = {"layer_min": 1, "layer_max": 1, "bottom": bottom, "top": top} expected_dims = {"index": 4, "species": 2} - return clip_arguments, expected_dims, does_not_raise() + return clip_arguments, expected_dims, does_not_raise(), does_not_raise() @staticmethod def case_clip_top_is_scalar(parameterizable_basic_dis): @@ -383,7 +410,7 @@ def case_clip_top_is_scalar(parameterizable_basic_dis): clip_arguments = {"layer_max": 2, "bottom": bottom, "top": top} expected_dims = {"index": 4, "species": 2} - return clip_arguments, expected_dims, does_not_raise() + return clip_arguments, expected_dims, does_not_raise(), does_not_raise() @staticmethod def case_clip_top_is_non_layered_structuredgrid(parameterizable_basic_dis): @@ -394,7 +421,7 @@ def case_clip_top_is_non_layered_structuredgrid(parameterizable_basic_dis): clip_arguments = {"layer_max": 2, "bottom": bottom, "top": top} expected_dims = {"index": 4, "species": 2} - return clip_arguments, expected_dims, does_not_raise() + return clip_arguments, expected_dims, does_not_raise(), does_not_raise() @staticmethod def case_clip_top_is_layered_structuredgrid(parameterizable_basic_dis): @@ -404,7 +431,7 @@ def case_clip_top_is_layered_structuredgrid(parameterizable_basic_dis): clip_arguments = {"layer_max": 2, "bottom": bottom, "top": top} expected_dims = {"index": 4, "species": 2} - return clip_arguments, expected_dims, does_not_raise() + return clip_arguments, expected_dims, does_not_raise(), does_not_raise() @staticmethod def case_clip_top_is_non_layered_unstructuredgrid(parameterizable_basic_dis): @@ -418,7 +445,7 @@ def case_clip_top_is_non_layered_unstructuredgrid(parameterizable_basic_dis): clip_arguments = {"layer_max": 2, "bottom": bottom, "top": top} expected_dims = {"index": 4, "species": 2} - return clip_arguments, expected_dims, does_not_raise() + return clip_arguments, expected_dims, does_not_raise(), does_not_raise() @staticmethod def case_clip_top_is_layered_unstructuredgrid(parameterizable_basic_dis): @@ -430,23 +457,23 @@ def case_clip_top_is_layered_unstructuredgrid(parameterizable_basic_dis): clip_arguments = {"layer_max": 2, "bottom": bottom, "top": top} expected_dims = {"index": 4, "species": 2} - return clip_arguments, expected_dims, does_not_raise() + return clip_arguments, expected_dims, does_not_raise(), does_not_raise() @staticmethod def case_clip_missing_top(parameterizable_basic_dis): _, _, bottom = parameterizable_basic_dis clip_arguments = {"layer_max": 2, "bottom": bottom} - expected_dims = {} - return clip_arguments, expected_dims, pytest.raises(ValueError) + expected_dims = {"index": 4, "species": 2} + return clip_arguments, expected_dims, pytest.raises(ValueError), does_not_raise() @staticmethod def case_clip_missing_bottom(parameterizable_basic_dis): _, top, _ = parameterizable_basic_dis clip_arguments = {"layer_max": 2, "top": top} - expected_dims = {} - return clip_arguments, expected_dims, pytest.raises(ValueError) + expected_dims = {"index": 4, "species": 2} + return clip_arguments, expected_dims, pytest.raises(ValueError), does_not_raise() @pytest.mark.parametrize( @@ -466,10 +493,10 @@ def case_clip_missing_bottom(parameterizable_basic_dis): indirect=True, ) @parametrize_with_cases( - ("clip_box_args", "expected_dims", "expectation"), cases=ClipBoxCases + ("clip_box_args", "expected_dims", "expectation", "_"), cases=ClipBoxCases ) -def test_clip_box__high_lvl_stationary( - well_high_lvl_test_data_stationary, clip_box_args, expected_dims, expectation +def test_clip_box__well_stationary( + well_high_lvl_test_data_stationary, clip_box_args, expected_dims, expectation, _ ): # Arrange wel = imod.mf6.Well(*well_high_lvl_test_data_stationary) @@ -482,6 +509,41 @@ def test_clip_box__high_lvl_stationary( assert dict(ds.dims) == expected_dims +@pytest.mark.parametrize( + "parameterizable_basic_dis", + [ + BasicDisSettings( + nlay=10, + zstop=-10.0, + xstart=50.0, + xstop=100.0, + ystart=50.0, + ystop=100.0, + nrow=10, + ncol=10, + ) + ], + indirect=True, +) +@parametrize_with_cases( + ("clip_box_args", "expected_dims", "_", "expectation"), cases=ClipBoxCases +) +def test_clip_box__layered_well_stationary( + well_high_lvl_test_data_stationary, clip_box_args, expected_dims, _, expectation +): + + x, y, _, _, rate_wel, concentration = well_high_lvl_test_data_stationary + layer = [1, 1, 1, 1, 9, 9, 9, 9] + wel = imod.mf6.LayeredWell(x, y, layer, rate_wel, concentration) + + with expectation: + # Act + ds = wel.clip_box(**clip_box_args).dataset + + # Assert + assert dict(ds.dims) == expected_dims + + @pytest.mark.parametrize( "parameterizable_basic_dis", [BasicDisSettings(nlay=10, zstop=-10.0)], @@ -550,7 +612,7 @@ def test_derive_cellid_from_points(basic_dis, well_high_lvl_test_data_stationary ) # Act - cellid = imod.mf6.wel.Well._Well__derive_cellid_from_points(idomain, x, y, layer) + cellid = imod.mf6.wel.Well._derive_cellid_from_points(idomain, x, y, layer) # Assert np.testing.assert_array_equal(cellid, cellid_expected) From aec7051ce0b7dc206ddc705b819a0c5d9737bea2 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Fri, 26 Jul 2024 16:22:26 +0200 Subject: [PATCH 07/25] Better naming for method --- imod/mf6/wel.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py index 7c0fd9e82..3870b3538 100644 --- a/imod/mf6/wel.py +++ b/imod/mf6/wel.py @@ -280,7 +280,7 @@ def to_mf6_pkg( wells_df = self._create_wells_df() nwells_df = len(wells_df["id"].unique()) - wells_assigned = self._create_assigned_wells(wells_df, active, top, bottom, k) + wells_assigned = self._assign_wells_to_layers(wells_df, active, top, bottom, k) nwells_assigned = ( 0 if wells_assigned.empty else len(wells_assigned["id"].unique()) @@ -336,7 +336,7 @@ def to_mf6_package_information(self, filtered_wells): def _create_wells_df(self): raise NotImplementedError("Method in abstract base class called") - def _create_assigned_wells( + def _assign_wells_to_layers( self, wells_df: pd.DataFrame, active: GridDataArray, @@ -630,7 +630,7 @@ def _validate(self, schemata: dict, **kwargs) -> dict[str, list[ValidationError] kwargs["screen_top"] = self.dataset["screen_top"] return Package._validate(self, schemata, **kwargs) - def _create_assigned_wells( + def _assign_wells_to_layers( self, wells_df: pd.DataFrame, active: GridDataArray, @@ -956,7 +956,7 @@ def clip_box( def _create_wells_df(self) -> pd.DataFrame: return self.dataset.to_dataframe() - def _create_assigned_wells( + def _assign_wells_to_layers( self, wells_df: pd.DataFrame, active: GridDataArray, From cb1ff38d306560101da303fd6e47bb52e9281c3c Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Fri, 26 Jul 2024 16:51:43 +0200 Subject: [PATCH 08/25] Update changelog --- docs/api/changelog.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst index 852a80f42..8bdd35b79 100644 --- a/docs/api/changelog.rst +++ b/docs/api/changelog.rst @@ -40,6 +40,8 @@ Added :class:`imod.mf6.HorizontalFlowBarrierResistance`, :class:`imod.mf6.HorizontalFlowBarrierMultiplier`, :class:`imod.mf6.HorizontalFlowBarrierHydraulicCharacteristic`. +- :class:`imod.mf6.LayeredWell` to specify wells directly to layers instead + assigning them with filter depths. [Unreleased] From 52c6c144004e88c6b5db719cb19f2e11736e5d92 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Fri, 26 Jul 2024 17:16:39 +0200 Subject: [PATCH 09/25] Fix example in docstring --- imod/mf6/wel.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py index 3870b3538..dcd8e77f8 100644 --- a/imod/mf6/wel.py +++ b/imod/mf6/wel.py @@ -816,7 +816,7 @@ class LayeredWell(GridAgnosticWell): >>> x = [81.0, 82.0] >>> rate = [1.0, 1.0] - >>> imod.mf6.Well(x, y, layer, rate) + >>> imod.mf6.LayeredWell(x, y, layer, rate) For a transient well: @@ -825,7 +825,7 @@ class LayeredWell(GridAgnosticWell): >>> rate_factor_time = xr.DataArray([0.5, 1.0], coords={"time": weltimes}, dims=("time",)) >>> rate_transient = rate_factor_time * xr.DataArray(rate, dims=("index",)) - >>> imod.mf6.Well(x, y, layer, rate_transient) + >>> imod.mf6.LayeredWell(x, y, layer, rate_transient) """ _pkg_id = "wel" From a66f54022ac97c00846ca95909f8f7955eb07709 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Fri, 26 Jul 2024 17:30:12 +0200 Subject: [PATCH 10/25] format --- imod/tests/test_mf6/test_mf6_wel.py | 34 +++++++++++++++++++---------- 1 file changed, 22 insertions(+), 12 deletions(-) diff --git a/imod/tests/test_mf6/test_mf6_wel.py b/imod/tests/test_mf6/test_mf6_wel.py index 85f1057cc..2dc0e0f3a 100644 --- a/imod/tests/test_mf6/test_mf6_wel.py +++ b/imod/tests/test_mf6/test_mf6_wel.py @@ -59,14 +59,16 @@ def case_well_stationary(self, well_high_lvl_test_data_stationary): return obj, dims_expected, cellid_expected, rate_expected def case_well_stationary_multilevel(self, well_high_lvl_test_data_stationary): - x, y, screen_top, _, rate_wel, concentration = well_high_lvl_test_data_stationary + x, y, screen_top, _, rate_wel, concentration = ( + well_high_lvl_test_data_stationary + ) screen_bottom = [-20.0] * 8 obj = imod.mf6.Well(x, y, screen_top, screen_bottom, rate_wel, concentration) dims_expected = { - "ncellid": 12, - "nmax_cellid": 3, - "species": 2, - } + "ncellid": 12, + "nmax_cellid": 3, + "species": 2, + } cellid_expected = np.array( [ [1, 1, 9], @@ -164,11 +166,10 @@ def case_layered_well_transient(self, well_high_lvl_test_data_transient): @parametrize_with_cases( - ["wel", "dims_expected", "cellid_expected", "rate_expected"], cases=GridAgnosticWellCases + ["wel", "dims_expected", "cellid_expected", "rate_expected"], + cases=GridAgnosticWellCases, ) -def test_to_mf6_pkg( - basic_dis, wel, dims_expected, cellid_expected, rate_expected -): +def test_to_mf6_pkg(basic_dis, wel, dims_expected, cellid_expected, rate_expected): # Arrange idomain, top, bottom = basic_dis active = idomain == 1 @@ -465,7 +466,12 @@ def case_clip_missing_top(parameterizable_basic_dis): clip_arguments = {"layer_max": 2, "bottom": bottom} expected_dims = {"index": 4, "species": 2} - return clip_arguments, expected_dims, pytest.raises(ValueError), does_not_raise() + return ( + clip_arguments, + expected_dims, + pytest.raises(ValueError), + does_not_raise(), + ) @staticmethod def case_clip_missing_bottom(parameterizable_basic_dis): @@ -473,7 +479,12 @@ def case_clip_missing_bottom(parameterizable_basic_dis): clip_arguments = {"layer_max": 2, "top": top} expected_dims = {"index": 4, "species": 2} - return clip_arguments, expected_dims, pytest.raises(ValueError), does_not_raise() + return ( + clip_arguments, + expected_dims, + pytest.raises(ValueError), + does_not_raise(), + ) @pytest.mark.parametrize( @@ -531,7 +542,6 @@ def test_clip_box__well_stationary( def test_clip_box__layered_well_stationary( well_high_lvl_test_data_stationary, clip_box_args, expected_dims, _, expectation ): - x, y, _, _, rate_wel, concentration = well_high_lvl_test_data_stationary layer = [1, 1, 1, 1, 9, 9, 9, 9] wel = imod.mf6.LayeredWell(x, y, layer, rate_wel, concentration) From b42fcd6b58d9fb6aba2cc87b4ea1d273902b2292 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Fri, 26 Jul 2024 18:29:31 +0200 Subject: [PATCH 11/25] Move grouping logic to prepare well rates to separate helper function. --- imod/mf6/wel.py | 117 +++++++++++++++++++++++++++++++----------------- 1 file changed, 75 insertions(+), 42 deletions(-) diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py index dcd8e77f8..90dc7c5bc 100644 --- a/imod/mf6/wel.py +++ b/imod/mf6/wel.py @@ -71,6 +71,30 @@ def mask_2D(package: GridAgnosticWell, domain_2d: GridDataArray) -> GridAgnostic return cls._from_dataset(selection) +def _prepare_well_rates_from_groups( + df_groups: pd.api.typing.DataFrameGroupBy, times: list[datetime] + ) -> xr.DataArray: + """ + Prepare well rates from dataframe groups, grouped by unique well locations. + """ + # Resample times per group + df_resampled_groups = [ + resample_timeseries(df_group, times) for df_group in df_groups + ] + # Convert dataframes all groups to DataArrays + da_groups = [ + xr.DataArray(df_group["rate"], dims=("time"), coords={"time": df_group["time"]}) + for df_group in df_resampled_groups + ] + # Assign index coordinates + da_groups = [ + da_group.expand_dims(dim="index").assign_coords(index=[i]) + for i, da_group in enumerate(da_groups) + ] + # Concatenate datarrays along index dimension + return xr.concat(da_groups, dim="index") + + class GridAgnosticWell(BoundaryCondition, IPointDataPackage, abc.ABC): """ Abstract base class for grid agnostic wells @@ -673,23 +697,20 @@ def from_imod5_data( minimum_k: float = 0.1, minimum_thickness: float = 1.0, ) -> "Well": - if "layer" in imod5_data[key].keys(): - if imod5_data[key]["layer"] != 0: - log_msg = textwrap.dedent( - f""" - In well {key} a layer was assigned, but this is not - supported. Assignment will be done based on filter_top and - filter_bottom, and the chosen layer ({imod5_data[key]["layer"]}) - will be ignored.""" - ) - logger.log( - loglevel=LogLevel.WARNING, message=log_msg, additional_depth=2 - ) - - if ( - "filt_top" not in imod5_data[key]["dataframe"].columns - or "filt_bot" not in imod5_data[key]["dataframe"].columns - ): + pkg_data = imod5_data[key] + if "layer" in pkg_data.keys() and (pkg_data["layer"] != 0): + log_msg = textwrap.dedent( + f""" + In well {key} a layer was assigned, but this is not + supported. Assignment will be done based on filter_top and + filter_bottom, and the chosen layer ({pkg_data["layer"]}) + will be ignored.""" + ) + logger.log(loglevel=LogLevel.WARNING, message=log_msg, additional_depth=2) + + df: pd.DataFrame = pkg_data["dataframe"] + + if "filt_top" not in df.columns or "filt_bot" not in df.columns: log_msg = textwrap.dedent( f""" In well {key} the filt_top and filt_bot columns were not both found; @@ -698,37 +719,13 @@ def from_imod5_data( logger.log(loglevel=LogLevel.ERROR, message=log_msg, additional_depth=2) raise ValueError(log_msg) - df: pd.DataFrame = imod5_data[key]["dataframe"] - # Groupby unique wells, to get dataframes per time. colnames_group = ["x", "y", "filt_top", "filt_bot", "id"] wel_index, df_groups = zip(*df.groupby(colnames_group)) - # resample per group - # Unpack wel indices by zipping x, y, filt_top, filt_bot, id = zip(*wel_index) - - # resample times per group - df_resampled_groups = [] - for df_group in df_groups: - df_group = resample_timeseries(df_group, times) - df_resampled_groups.append(df_group) - - # Convert dataframes all groups to DataArrays - da_groups = [ - xr.DataArray( - df_group["rate"], dims=("time"), coords={"time": df_group["time"]} - ) - for df_group in df_resampled_groups - ] - # Assign index coordinates - da_groups = [ - da_group.expand_dims(dim="index").assign_coords(index=[i]) - for i, da_group in enumerate(da_groups) - ] - # Concatenate datarrays along index dimension - well_rate = xr.concat(da_groups, dim="index") + well_rate = _prepare_well_rates_from_groups(df_groups, times) return cls( x=np.array(x, dtype=float), @@ -966,6 +963,42 @@ def _assign_wells_to_layers( ): return wells_df + @classmethod + def from_imod5_data( + cls, + key: str, + imod5_data: dict[str, dict[str, GridDataArray]], + times: list[datetime], + minimum_k: float = 0.1, + minimum_thickness: float = 1.0, + ) -> "LayeredWell": + pkg_data = imod5_data[key] + + if ("layer" not in pkg_data.keys()) or (pkg_data["layer"] == 0): + log_msg = textwrap.dedent( + f"""In well {key} no layer was assigned, but this is required.""" + ) + logger.log(loglevel=LogLevel.ERROR, message=log_msg, additional_depth=2) + + df: pd.DataFrame = pkg_data["dataframe"] + + # Groupby unique wells, to get dataframes per time. + colnames_group = ["x", "y", "layer", "id"] + wel_index, df_groups = zip(*df.groupby(colnames_group)) + + # Unpack wel indices by zipping + x, y, layer, id = zip(*wel_index) + well_rate = _prepare_well_rates_from_groups(df_groups, times) + + return cls( + x=np.array(x, dtype=float), + y=np.array(y, dtype=float), + layer=np.array(layer, dtype=int), + rate=well_rate, + minimum_k=minimum_k, + minimum_thickness=minimum_thickness, + ) + class WellDisStructured(DisStructuredBoundaryCondition): """ From 846faaa3bd5e1dc97cdf16c1f14a0585695fb166 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Fri, 26 Jul 2024 18:41:34 +0200 Subject: [PATCH 12/25] If layer not 0, make LayeredWell --- imod/mf6/model_gwf.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/imod/mf6/model_gwf.py b/imod/mf6/model_gwf.py index b479827cf..f94340a89 100644 --- a/imod/mf6/model_gwf.py +++ b/imod/mf6/model_gwf.py @@ -31,7 +31,7 @@ from imod.mf6.riv import River from imod.mf6.sto import StorageCoefficient from imod.mf6.utilities.regrid import RegridderWeightsCache -from imod.mf6.wel import Well +from imod.mf6.wel import LayeredWell, Well from imod.prepare.topsystem.default_allocation_methods import ( SimulationAllocationOptions, SimulationDistributingOptions, @@ -281,7 +281,11 @@ def from_imod5_data( imod5_keys = list(imod5_data.keys()) wel_keys = [key for key in imod5_keys if key[0:3] == "wel"] for wel_key in wel_keys: - result[wel_key] = Well.from_imod5_data(wel_key, imod5_data, times) + layer = imod5_data[wel_key]["layer"] + if layer == 0: + result[wel_key] = Well.from_imod5_data(wel_key, imod5_data, times) + else: + result[wel_key] = LayeredWell.from_imod5_data(wel_key, imod5_data, times) # import drainage From 46309f111bd1959c020bee0c5c7283bcfe8fab21 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Fri, 26 Jul 2024 18:43:09 +0200 Subject: [PATCH 13/25] format --- imod/mf6/model_gwf.py | 4 +++- imod/mf6/wel.py | 4 ++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/imod/mf6/model_gwf.py b/imod/mf6/model_gwf.py index f94340a89..ba154b483 100644 --- a/imod/mf6/model_gwf.py +++ b/imod/mf6/model_gwf.py @@ -285,7 +285,9 @@ def from_imod5_data( if layer == 0: result[wel_key] = Well.from_imod5_data(wel_key, imod5_data, times) else: - result[wel_key] = LayeredWell.from_imod5_data(wel_key, imod5_data, times) + result[wel_key] = LayeredWell.from_imod5_data( + wel_key, imod5_data, times + ) # import drainage diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py index 90dc7c5bc..86a2146cb 100644 --- a/imod/mf6/wel.py +++ b/imod/mf6/wel.py @@ -72,8 +72,8 @@ def mask_2D(package: GridAgnosticWell, domain_2d: GridDataArray) -> GridAgnostic def _prepare_well_rates_from_groups( - df_groups: pd.api.typing.DataFrameGroupBy, times: list[datetime] - ) -> xr.DataArray: + df_groups: pd.api.typing.DataFrameGroupBy, times: list[datetime] +) -> xr.DataArray: """ Prepare well rates from dataframe groups, grouped by unique well locations. """ From 88158f74b5758ff925e4e9b6f839c762725fc632 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Mon, 29 Jul 2024 10:38:53 +0200 Subject: [PATCH 14/25] Better variable naming --- imod/mf6/wel.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py index 86a2146cb..6ba78d100 100644 --- a/imod/mf6/wel.py +++ b/imod/mf6/wel.py @@ -72,14 +72,14 @@ def mask_2D(package: GridAgnosticWell, domain_2d: GridDataArray) -> GridAgnostic def _prepare_well_rates_from_groups( - df_groups: pd.api.typing.DataFrameGroupBy, times: list[datetime] + unique_well_groups: pd.api.typing.DataFrameGroupBy, times: list[datetime] ) -> xr.DataArray: """ Prepare well rates from dataframe groups, grouped by unique well locations. """ # Resample times per group df_resampled_groups = [ - resample_timeseries(df_group, times) for df_group in df_groups + resample_timeseries(df_group, times) for df_group in unique_well_groups ] # Convert dataframes all groups to DataArrays da_groups = [ @@ -721,11 +721,11 @@ def from_imod5_data( # Groupby unique wells, to get dataframes per time. colnames_group = ["x", "y", "filt_top", "filt_bot", "id"] - wel_index, df_groups = zip(*df.groupby(colnames_group)) + wel_index, unique_well_groups = zip(*df.groupby(colnames_group)) # Unpack wel indices by zipping x, y, filt_top, filt_bot, id = zip(*wel_index) - well_rate = _prepare_well_rates_from_groups(df_groups, times) + well_rate = _prepare_well_rates_from_groups(unique_well_groups, times) return cls( x=np.array(x, dtype=float), @@ -984,11 +984,11 @@ def from_imod5_data( # Groupby unique wells, to get dataframes per time. colnames_group = ["x", "y", "layer", "id"] - wel_index, df_groups = zip(*df.groupby(colnames_group)) + wel_index, unique_well_groups = zip(*df.groupby(colnames_group)) # Unpack wel indices by zipping x, y, layer, id = zip(*wel_index) - well_rate = _prepare_well_rates_from_groups(df_groups, times) + well_rate = _prepare_well_rates_from_groups(unique_well_groups, times) return cls( x=np.array(x, dtype=float), From efd930cc0ae52d151af437311befdb4b4dc8f3c6 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Mon, 29 Jul 2024 14:22:48 +0200 Subject: [PATCH 15/25] isinstance GridAgnosticWell instead of Well --- imod/mf6/model.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/imod/mf6/model.py b/imod/mf6/model.py index 031c7623d..a30b6291b 100644 --- a/imod/mf6/model.py +++ b/imod/mf6/model.py @@ -27,6 +27,7 @@ from imod.mf6.utilities.mf6hfb import merge_hfb_packages from imod.mf6.utilities.regrid import RegridderWeightsCache, _regrid_like from imod.mf6.validation import pkg_errors_to_status_info +from imod.mf6.wel import GridAgnosticWell from imod.mf6.write_context import WriteContext from imod.schemata import ValidationError from imod.typing import GridDataArray @@ -267,7 +268,7 @@ def write( mf6_hfb_ls: List[HorizontalFlowBarrierBase] = [] for pkg_name, pkg in self.items(): try: - if isinstance(pkg, imod.mf6.Well): + if isinstance(pkg, GridAgnosticWell): top, bottom, idomain = self.__get_domain_geometry() k = self.__get_k() mf6_well_pkg = pkg.to_mf6_pkg( From e40f91e249145b36bd36c3357877aeb6e3ed3022 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Mon, 29 Jul 2024 14:23:10 +0200 Subject: [PATCH 16/25] Add layer to dataframe --- imod/mf6/wel.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py index 6ba78d100..4505c417c 100644 --- a/imod/mf6/wel.py +++ b/imod/mf6/wel.py @@ -982,6 +982,9 @@ def from_imod5_data( df: pd.DataFrame = pkg_data["dataframe"] + # Add layer to dataframe. + df["layer"] = pkg_data["layer"] + # Groupby unique wells, to get dataframes per time. colnames_group = ["x", "y", "layer", "id"] wel_index, unique_well_groups = zip(*df.groupby(colnames_group)) From 58286abc9a3c41b2a41a3edbb34bd908edd90698 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Mon, 29 Jul 2024 14:25:09 +0200 Subject: [PATCH 17/25] Parametrize test with LayeredWell case as well. --- imod/tests/test_mf6/test_mf6_wel.py | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/imod/tests/test_mf6/test_mf6_wel.py b/imod/tests/test_mf6/test_mf6_wel.py index 2dc0e0f3a..dfa7304e7 100644 --- a/imod/tests/test_mf6/test_mf6_wel.py +++ b/imod/tests/test_mf6/test_mf6_wel.py @@ -10,7 +10,7 @@ import pytest import xarray as xr import xugrid as xu -from pytest_cases import parametrize_with_cases +from pytest_cases import parametrize, parametrize_with_cases import imod from imod.formats.prj.prj import open_projectfile_data @@ -19,7 +19,7 @@ from imod.mf6.dis import StructuredDiscretization from imod.mf6.npf import NodePropertyFlow from imod.mf6.utilities.grid import broadcast_to_full_domain -from imod.mf6.wel import Well +from imod.mf6.wel import LayeredWell, Well from imod.mf6.write_context import WriteContext from imod.schemata import ValidationError from imod.tests.fixtures.flow_basic_fixture import BasicDisSettings @@ -901,8 +901,9 @@ def test_render__concentration_dis_vertices_transient(well_test_data_transient): ) # check salinity and temperature was written to period data +@parametrize("wel_class", [Well, LayeredWell]) @pytest.mark.usefixtures("imod5_dataset") -def test_import_and_convert_to_mf6(imod5_dataset, tmp_path): +def test_import_and_convert_to_mf6(imod5_dataset, tmp_path, wel_class): data = imod5_dataset[0] target_dis = StructuredDiscretization.from_imod5_data(data) target_npf = NodePropertyFlow.from_imod5_data(data, target_dis.dataset["idomain"]) @@ -910,7 +911,7 @@ def test_import_and_convert_to_mf6(imod5_dataset, tmp_path): times = list(pd.date_range(datetime(1989, 1, 1), datetime(2013, 1, 1), 8400)) # import grid-agnostic well from imod5 data (it contains 1 well) - wel = Well.from_imod5_data("wel-1", data, times) + wel = wel_class.from_imod5_data("wel-1", data, times) assert wel.dataset["x"].values[0] == 197910.0 assert wel.dataset["y"].values[0] == 362860.0 assert np.mean(wel.dataset["rate"].values) == -317.1681465863781 @@ -932,8 +933,9 @@ def test_import_and_convert_to_mf6(imod5_dataset, tmp_path): mf6_well.write("wel", [], write_context) +@parametrize("wel_class", [Well, LayeredWell]) @pytest.mark.usefixtures("well_regular_import_prj") -def test_import_multiple_wells(well_regular_import_prj): +def test_import_multiple_wells(well_regular_import_prj, wel_class): imod5dict = open_projectfile_data(well_regular_import_prj) times = [ datetime(1981, 11, 30), @@ -945,8 +947,8 @@ def test_import_multiple_wells(well_regular_import_prj): ] # import grid-agnostic well from imod5 data (it contains 2 packages with 3 wells each) - wel1 = imod.mf6.Well.from_imod5_data("wel-1", imod5dict[0], times) - wel2 = imod.mf6.Well.from_imod5_data("wel-2", imod5dict[0], times) + wel1 = wel_class.from_imod5_data("wel-1", imod5dict[0], times) + wel2 = wel_class.from_imod5_data("wel-2", imod5dict[0], times) assert np.all(wel1.x == np.array([191112.11, 191171.96, 191231.52])) assert np.all(wel2.x == np.array([191112.11, 191171.96, 191231.52])) @@ -954,8 +956,9 @@ def test_import_multiple_wells(well_regular_import_prj): assert wel2.dataset["rate"].shape == (6, 3) +@parametrize("wel_class", [Well, LayeredWell]) @pytest.mark.usefixtures("well_duplication_import_prj") -def test_import_from_imod5_with_duplication(well_duplication_import_prj): +def test_import_from_imod5_with_duplication(well_duplication_import_prj, wel_class): imod5dict = open_projectfile_data(well_duplication_import_prj) times = [ datetime(1981, 11, 30), @@ -966,8 +969,8 @@ def test_import_from_imod5_with_duplication(well_duplication_import_prj): datetime(1982, 4, 30), ] # import grid-agnostic well from imod5 data (it contains 2 packages with 3 wells each) - wel1 = imod.mf6.Well.from_imod5_data("wel-1", imod5dict[0], times) - wel2 = imod.mf6.Well.from_imod5_data("wel-2", imod5dict[0], times) + wel1 = wel_class.from_imod5_data("wel-1", imod5dict[0], times) + wel2 = wel_class.from_imod5_data("wel-2", imod5dict[0], times) assert np.all(wel1.x == np.array([191171.96, 191231.52, 191231.52])) assert np.all(wel2.x == np.array([191112.11, 191171.96, 191231.52])) From b5b261c68361050d61e32f3639a41c9a4e154118 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Mon, 29 Jul 2024 14:25:34 +0200 Subject: [PATCH 18/25] Add test to verify the right classes are instantiated. --- imod/tests/test_mf6/test_mf6_simulation.py | 33 +++++++++++++++++++++- 1 file changed, 32 insertions(+), 1 deletion(-) diff --git a/imod/tests/test_mf6/test_mf6_simulation.py b/imod/tests/test_mf6/test_mf6_simulation.py index 88d123169..893ffe9bf 100644 --- a/imod/tests/test_mf6/test_mf6_simulation.py +++ b/imod/tests/test_mf6/test_mf6_simulation.py @@ -17,6 +17,7 @@ import xugrid as xu import imod +from imod.mf6 import LayeredWell, Well from imod.mf6.model import Modflow6Model from imod.mf6.multimodel.modelsplitter import PartitionInfo from imod.mf6.oc import OutputControl @@ -507,7 +508,37 @@ def test_import_from_imod5(imod5_dataset, tmp_path): @pytest.mark.usefixtures("imod5_dataset") -def test_import_from_imod5_nonstandard_regridding(imod5_dataset, tmp_path): +def test_import_from_imod5__correct_well_type(imod5_dataset): + # Unpack + imod5_data = imod5_dataset[0] + period_data = imod5_dataset[1] + # Temporarily change layer number to 0, to force Well object instead of + # LayeredWell + original_wel_layer = imod5_data["wel-1"]["layer"] + imod5_data["wel-1"]["layer"] = 0 + # Other arrangement + default_simulation_allocation_options = SimulationAllocationOptions + default_simulation_distributing_options = SimulationDistributingOptions + datelist = pd.date_range(start="1/1/1989", end="1/1/2013", freq="W") + + # Act + simulation = Modflow6Simulation.from_imod5_data( + imod5_data, + period_data, + default_simulation_allocation_options, + default_simulation_distributing_options, + datelist, + ) + # Set layer back to right value (before AssertionError might be thrown) + imod5_data["wel-1"]["layer"] = original_wel_layer + # Assert + assert isinstance(simulation["imported_model"]["wel-1"], Well) + assert isinstance(simulation["imported_model"]["wel-2"], LayeredWell) + assert isinstance(simulation["imported_model"]["wel-3"], LayeredWell) + + +@pytest.mark.usefixtures("imod5_dataset") +def test_import_from_imod5__nonstandard_regridding(imod5_dataset, tmp_path): imod5_data = imod5_dataset[0] period_data = imod5_dataset[1] default_simulation_allocation_options = SimulationAllocationOptions From 30bd4ea13b568a47960b4196b41e633e2cd50b03 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Mon, 29 Jul 2024 17:19:33 +0200 Subject: [PATCH 19/25] Test if subclass explicitly --- imod/mf6/model.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/imod/mf6/model.py b/imod/mf6/model.py index a30b6291b..d3bf56cee 100644 --- a/imod/mf6/model.py +++ b/imod/mf6/model.py @@ -268,7 +268,7 @@ def write( mf6_hfb_ls: List[HorizontalFlowBarrierBase] = [] for pkg_name, pkg in self.items(): try: - if isinstance(pkg, GridAgnosticWell): + if issubclass(type(pkg), GridAgnosticWell): top, bottom, idomain = self.__get_domain_geometry() k = self.__get_k() mf6_well_pkg = pkg.to_mf6_pkg( @@ -284,7 +284,7 @@ def write( globaltimes=globaltimes, write_context=pkg_write_context, ) - elif isinstance(pkg, imod.mf6.HorizontalFlowBarrierBase): + elif issubclass(type(pkg), imod.mf6.HorizontalFlowBarrierBase): mf6_hfb_ls.append(pkg) else: pkg.write( From 88f934c4c0fcc6d6afe8ed318ae9cfeea6d53d48 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Mon, 29 Jul 2024 17:20:32 +0200 Subject: [PATCH 20/25] Remove code reduncancy where index was retrieved from other df. --- imod/mf6/wel.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py index 4505c417c..82d185582 100644 --- a/imod/mf6/wel.py +++ b/imod/mf6/wel.py @@ -123,7 +123,7 @@ def _create_cellid(self, wells_assigned: pd.DataFrame, active: xr.DataArray): return self._derive_cellid_from_points(like, **d_for_cellid) def _create_dataset_vars( - self, wells_assigned: pd.DataFrame, wells_df: pd.DataFrame, cellid: xr.DataArray + self, wells_assigned: pd.DataFrame, cellid: xr.DataArray ) -> xr.Dataset: """ Create dataset with all variables (rate, concentration), with a similar shape as the cellids. @@ -136,7 +136,7 @@ def _create_dataset_vars( # "rate" variable in conversion from multi-indexed DataFrame to xarray # DataArray results in duplicated values for "rate" along dimension # "species". Select first species to reduce this again. - index_names = wells_df.index.names + index_names = wells_assigned.index.names if "species" in index_names: ds_vars["rate"] = ds_vars["rate"].isel(species=0) @@ -321,7 +321,7 @@ def to_mf6_pkg( ds = xr.Dataset() ds["cellid"] = self._create_cellid(wells_assigned, active) - ds_vars = self._create_dataset_vars(wells_assigned, wells_df, ds["cellid"]) + ds_vars = self._create_dataset_vars(wells_assigned, ds["cellid"]) ds = ds.assign(**ds_vars.data_vars) ds = remove_inactive(ds, active) From d5ff265f745ce4e16a076f9a119bcc356c80988a Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Mon, 29 Jul 2024 17:27:13 +0200 Subject: [PATCH 21/25] Specify return type --- imod/mf6/wel.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py index 82d185582..53846b03f 100644 --- a/imod/mf6/wel.py +++ b/imod/mf6/wel.py @@ -367,7 +367,7 @@ def _assign_wells_to_layers( top: GridDataArray, bottom: GridDataArray, k: GridDataArray, - ): + ) -> pd.DataFrame: raise NotImplementedError("Method in abstract base class called") @@ -661,7 +661,7 @@ def _assign_wells_to_layers( top: GridDataArray, bottom: GridDataArray, k: GridDataArray, - ): + ) -> pd.DataFrame: # Ensure top, bottom & k # are broadcasted to 3d grid like = ones_like(active) @@ -960,7 +960,7 @@ def _assign_wells_to_layers( top: GridDataArray, bottom: GridDataArray, k: GridDataArray, - ): + ) -> pd.DataFrame: return wells_df @classmethod From 42c85d7e73644a594ebb9ba8045f69e7ecc5cc55 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Mon, 29 Jul 2024 17:33:52 +0200 Subject: [PATCH 22/25] Add return types --- imod/mf6/wel.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py index 53846b03f..8c7ecb97f 100644 --- a/imod/mf6/wel.py +++ b/imod/mf6/wel.py @@ -112,7 +112,9 @@ def y(self) -> npt.NDArray[np.float64]: def is_grid_agnostic_package(cls) -> bool: return True - def _create_cellid(self, wells_assigned: pd.DataFrame, active: xr.DataArray): + def _create_cellid( + self, wells_assigned: pd.DataFrame, active: xr.DataArray + ) -> GridDataArray: like = ones_like(active) # Groupby index and select first, to unset any duplicate records @@ -340,7 +342,7 @@ def to_mf6_pkg( return Mf6Wel(**ds.data_vars) - def to_mf6_package_information(self, filtered_wells): + def to_mf6_package_information(self, filtered_wells: pd.DataFrame) -> str: message = textwrap.dedent( """Some wells were not placed in the MF6 well package. This can be due to inactive cells or permeability/thickness constraints.\n""" @@ -357,7 +359,7 @@ def to_mf6_package_information(self, filtered_wells): message += f" id = {ids} x = {x} y = {y} \n" return message - def _create_wells_df(self): + def _create_wells_df(self) -> pd.DataFrame: raise NotImplementedError("Method in abstract base class called") def _assign_wells_to_layers( From 77f6a243edde68eee373a8467fc90074984135ff Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Mon, 29 Jul 2024 17:40:13 +0200 Subject: [PATCH 23/25] Expand log message with info that imod.mf6.LayeredWell should be used. --- imod/mf6/wel.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py index 8c7ecb97f..83919934a 100644 --- a/imod/mf6/wel.py +++ b/imod/mf6/wel.py @@ -703,10 +703,11 @@ def from_imod5_data( if "layer" in pkg_data.keys() and (pkg_data["layer"] != 0): log_msg = textwrap.dedent( f""" - In well {key} a layer was assigned, but this is not - supported. Assignment will be done based on filter_top and - filter_bottom, and the chosen layer ({pkg_data["layer"]}) - will be ignored.""" + In well {key} a layer was assigned, but this is not supported. + Assignment will be done based on filter_top and filter_bottom, + and the chosen layer ({pkg_data["layer"]}) will be ignored. To + specify by layer, use imod.mf6.LayeredWell. + """ ) logger.log(loglevel=LogLevel.WARNING, message=log_msg, additional_depth=2) @@ -715,8 +716,10 @@ def from_imod5_data( if "filt_top" not in df.columns or "filt_bot" not in df.columns: log_msg = textwrap.dedent( f""" - In well {key} the filt_top and filt_bot columns were not both found; - this is not supported for import.""" + In well {key} the filt_top and filt_bot columns were not both + found; this is not supported for import. To specify by layer, + use imod.mf6.LayeredWell. + """ ) logger.log(loglevel=LogLevel.ERROR, message=log_msg, additional_depth=2) raise ValueError(log_msg) From 16b4798924b3b27f2dc73313343c871c11a8314b Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Mon, 29 Jul 2024 17:43:26 +0200 Subject: [PATCH 24/25] Add LayeredWell class to package sanity --- imod/tests/test_mf6/test_package_sanity.py | 1 + 1 file changed, 1 insertion(+) diff --git a/imod/tests/test_mf6/test_package_sanity.py b/imod/tests/test_mf6/test_package_sanity.py index 545eb3d4a..cfd2c0851 100644 --- a/imod/tests/test_mf6/test_package_sanity.py +++ b/imod/tests/test_mf6/test_package_sanity.py @@ -42,6 +42,7 @@ imod.mf6.HorizontalFlowBarrierHydraulicCharacteristic, imod.mf6.HorizontalFlowBarrierMultiplier, imod.mf6.HorizontalFlowBarrierResistance, + imod.mf6.LayeredWell, ] From 5477bef84c9f6eee6beaef684c644fd6bc166fd2 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Mon, 29 Jul 2024 18:07:14 +0200 Subject: [PATCH 25/25] Update docstring, scalars not supported in docstring. --- imod/mf6/wel.py | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py index 83919934a..fb5195b0a 100644 --- a/imod/mf6/wel.py +++ b/imod/mf6/wel.py @@ -384,15 +384,15 @@ class Well(GridAgnosticWell): Parameters ---------- - y: float or list of floats or np.array of floats + y: list of floats or np.array of floats is the y location of the well. - x: float or list of floats or np.array of floats + x: list of floats or np.array of floats is the x location of the well. - screen_top: float or list of floats or np.array of floats + screen_top: list of floats or np.array of floats is the top of the well screen. - screen_bottom: float or list of floats or np.array of floats + screen_bottom: list of floats or np.array of floats is the bottom of the well screen. - rate: float, list of floats or xr.DataArray + rate: list of floats or xr.DataArray is the volumetric well rate. A positive value indicates well (injection) and a negative value indicates discharge (extraction) (q). If provided as DataArray, an ``"index"`` dimension is required and an @@ -747,20 +747,21 @@ class LayeredWell(GridAgnosticWell): """ Agnostic WEL package, which accepts x, y and layers. - This package can be written to any provided model grid. - Any number of WEL Packages can be specified for a single groundwater flow model. + This package can be written to any provided model grid, given that it has + enough layers. Any number of WEL Packages can be specified for a single + groundwater flow model. https://water.usgs.gov/water-resources/software/MODFLOW-6/mf6io_6.0.4.pdf#page=63 Parameters ---------- - y: float or list of floats or np.array of floats + y: list of floats or np.array of floats is the y location of the well. - x: float or list of floats or np.array of floats + x: list of floats or np.array of floats is the x location of the well. - layer: int or list of ints or np.array of ints + layer: list of ints or np.array of ints is the layer of the well. - rate: float, list of floats or xr.DataArray + rate: list of floats or xr.DataArray is the volumetric well rate. A positive value indicates well (injection) and a negative value indicates discharge (extraction) (q). If provided as DataArray, an ``"index"`` dimension is required and an