Coverage for C:\src\imod-python\imod\mf6\wel.py: 95%
227 statements
« prev ^ index » next coverage.py v7.4.4, created at 2024-04-08 13:27 +0200
« prev ^ index » next coverage.py v7.4.4, created at 2024-04-08 13:27 +0200
1from __future__ import annotations
3import warnings
4from typing import Any, Optional, Tuple, Union
6import cftime
7import numpy as np
8import numpy.typing as npt
9import pandas as pd
10import xarray as xr
11import xugrid as xu
13import imod
14from imod.logging import init_log_decorator
15from imod.mf6.boundary_condition import (
16 BoundaryCondition,
17 DisStructuredBoundaryCondition,
18 DisVerticesBoundaryCondition,
19)
20from imod.mf6.interfaces.ipointdatapackage import IPointDataPackage
21from imod.mf6.mf6_wel_adapter import Mf6Wel
22from imod.mf6.package import Package
23from imod.mf6.utilities.dataset import remove_inactive
24from imod.mf6.utilities.regrid import RegridderType
25from imod.mf6.validation import validation_pkg_error_message
26from imod.mf6.write_context import WriteContext
27from imod.prepare import assign_wells
28from imod.prepare.layer import create_layered_top
29from imod.schemata import (
30 AnyNoDataSchema,
31 DTypeSchema,
32 EmptyIndexesSchema,
33 ValidationError,
34)
35from imod.select.points import points_indices, points_values
36from imod.typing import GridDataArray
37from imod.typing.grid import is_spatial_2D, ones_like
38from imod.util.structured import values_within_range
41def _assign_dims(arg: Any) -> Tuple | xr.DataArray:
42 is_da = isinstance(arg, xr.DataArray)
43 if is_da and "time" in arg.coords:
44 if arg.ndim != 2:
45 raise ValueError("time varying variable: must be 2d")
46 if arg.dims[0] != "time":
47 arg = arg.transpose()
48 da = xr.DataArray(
49 data=arg.values, coords={"time": arg["time"]}, dims=["time", "index"]
50 )
51 return da
52 elif is_da:
53 return "index", arg.values
54 else:
55 return "index", arg
58def mask_2D(package: Well, domain_2d: GridDataArray) -> Well:
59 point_active = points_values(domain_2d, x=package.x, y=package.y)
61 is_inside_exterior = point_active == 1
62 selection = package.dataset.loc[{"index": is_inside_exterior}]
64 cls = type(package)
65 new = cls.__new__(cls)
66 new.dataset = selection
67 return new
70class Well(BoundaryCondition, IPointDataPackage):
71 """
72 Agnostic WEL package, which accepts x, y and a top and bottom of the well screens.
74 This package can be written to any provided model grid.
75 Any number of WEL Packages can be specified for a single groundwater flow model.
76 https://water.usgs.gov/water-resources/software/MODFLOW-6/mf6io_6.0.4.pdf#page=63
78 Parameters
79 ----------
81 y: float or list of floats
82 is the y location of the well.
83 x: float or list of floats
84 is the x location of the well.
85 screen_top: float or list of floats
86 is the top of the well screen.
87 screen_bottom: float or list of floats
88 is the bottom of the well screen.
89 rate: float, list of floats or xr.DataArray
90 is the volumetric well rate. A positive value indicates well
91 (injection) and a negative value indicates discharge (extraction) (q).
92 If provided as DataArray, an ``"index"`` dimension is required and an
93 optional ``"time"`` dimension and coordinate specify transient input.
94 In the latter case, it is important that dimensions are in the order:
95 ``("time", "index")``
96 concentration: array of floats (xr.DataArray, optional)
97 if this flow package is used in simulations also involving transport, then this array is used
98 as the concentration for inflow over this boundary.
99 concentration_boundary_type: ({"AUX", "AUXMIXED"}, optional)
100 if this flow package is used in simulations also involving transport, then this keyword specifies
101 how outflow over this boundary is computed.
102 id: list of Any, optional
103 assign an identifier code to each well. if not provided, one will be generated
104 Must be convertible to string, and unique entries.
105 minimum_k: float, optional
106 on creating point wells, no point wells will be placed in cells with a lower horizontal conductivity than this
107 minimum_thickness: float, optional
108 on creating point wells, no point wells will be placed in cells with a lower thickness than this
109 print_input: ({True, False}, optional)
110 keyword to indicate that the list of well information will be written to
111 the listing file immediately after it is read.
112 Default is False.
113 print_flows: ({True, False}, optional)
114 Indicates that the list of well flow rates will be printed to the
115 listing file for every stress period time step in which "BUDGET PRINT"
116 is specified in Output Control. If there is no Output Control option
117 and PRINT FLOWS is specified, then flow rates are printed for the last
118 time step of each stress period.
119 Default is False.
120 save_flows: ({True, False}, optional)
121 Indicates that well flow terms will be written to the file specified
122 with "BUDGET FILEOUT" in Output Control.
123 Default is False.
124 observations: [Not yet supported.]
125 Default is None.
126 validate: {True, False}
127 Flag to indicate whether the package should be validated upon
128 initialization. This raises a ValidationError if package input is
129 provided in the wrong manner. Defaults to True.
130 repeat_stress: Optional[xr.DataArray] of datetimes
131 Used to repeat data for e.g. repeating stress periods such as
132 seasonality without duplicating the values. The DataArray should have
133 dimensions ``("repeat", "repeat_items")``. The ``repeat_items``
134 dimension should have size 2: the first value is the "key", the second
135 value is the "value". For the "key" datetime, the data of the "value"
136 datetime will be used. Can also be set with a dictionary using the
137 ``set_repeat_stress`` method.
139 Examples
140 ---------
142 >>> screen_top = [0.0, 0.0]
143 >>> screen_bottom = [-2.0, -2.0]
144 >>> y = [83.0, 77.0]
145 >>> x = [81.0, 82.0]
146 >>> rate = [1.0, 1.0]
148 >>> imod.mf6.Well(x, y, screen_top, screen_bottom, rate)
150 For a transient well:
152 >>> weltimes = pd.date_range("2000-01-01", "2000-01-03")
154 >>> rate_factor_time = xr.DataArray([0.5, 1.0], coords={"time": weltimes}, dims=("time",))
155 >>> rate_transient = rate_factor_time * xr.DataArray(rate, dims=("index",))
157 >>> imod.mf6.Well(x, y, screen_top, screen_bottom, rate_transient)
158 """
160 @property
161 def x(self) -> npt.NDArray[np.float64]:
162 return self.dataset["x"].values
164 @property
165 def y(self) -> npt.NDArray[np.float64]:
166 return self.dataset["y"].values
168 _pkg_id = "wel"
170 _auxiliary_data = {"concentration": "species"}
171 _init_schemata = {
172 "screen_top": [DTypeSchema(np.floating)],
173 "screen_bottom": [DTypeSchema(np.floating)],
174 "y": [DTypeSchema(np.floating)],
175 "x": [DTypeSchema(np.floating)],
176 "rate": [DTypeSchema(np.floating)],
177 "concentration": [DTypeSchema(np.floating)],
178 }
179 _write_schemata = {
180 "screen_top": [AnyNoDataSchema(), EmptyIndexesSchema()],
181 "screen_bottom": [AnyNoDataSchema(), EmptyIndexesSchema()],
182 "y": [AnyNoDataSchema(), EmptyIndexesSchema()],
183 "x": [AnyNoDataSchema(), EmptyIndexesSchema()],
184 "rate": [AnyNoDataSchema(), EmptyIndexesSchema()],
185 "concentration": [AnyNoDataSchema(), EmptyIndexesSchema()],
186 }
188 _regrid_method: dict[str, Tuple[RegridderType, str]] = {}
190 @init_log_decorator()
191 def __init__(
192 self,
193 x: list[float],
194 y: list[float],
195 screen_top: list[float],
196 screen_bottom: list[float],
197 rate: list[float] | xr.DataArray,
198 concentration: Optional[list[float] | xr.DataArray] = None,
199 concentration_boundary_type="aux",
200 id: Optional[list[Any]] = None,
201 minimum_k: float = 0.1,
202 minimum_thickness: float = 1.0,
203 print_input: bool = False,
204 print_flows: bool = False,
205 save_flows: bool = False,
206 observations=None,
207 validate: bool = True,
208 repeat_stress: Optional[xr.DataArray] = None,
209 ):
210 if id is None:
211 id = [str(i) for i in range(len(x))]
212 else:
213 set_id = set(id)
214 if len(id) != len(set_id):
215 raise ValueError("id's must be unique")
216 id = [str(i) for i in id]
217 dict_dataset = {
218 "screen_top": _assign_dims(screen_top),
219 "screen_bottom": _assign_dims(screen_bottom),
220 "y": _assign_dims(y),
221 "x": _assign_dims(x),
222 "rate": _assign_dims(rate),
223 "id": _assign_dims(id),
224 "minimum_k": minimum_k,
225 "minimum_thickness": minimum_thickness,
226 "print_input": print_input,
227 "print_flows": print_flows,
228 "save_flows": save_flows,
229 "observations": observations,
230 "repeat_stress": repeat_stress,
231 "concentration": concentration,
232 "concentration_boundary_type": concentration_boundary_type,
233 }
234 super().__init__(dict_dataset)
235 # Set index as coordinate
236 index_coord = np.arange(self.dataset.dims["index"])
237 self.dataset = self.dataset.assign_coords(index=index_coord)
238 self._validate_init_schemata(validate)
240 @classmethod
241 def is_grid_agnostic_package(cls) -> bool:
242 return True
244 def clip_box(
245 self,
246 time_min: Optional[cftime.datetime | np.datetime64 | str] = None,
247 time_max: Optional[cftime.datetime | np.datetime64 | str] = None,
248 layer_min: Optional[int] = None,
249 layer_max: Optional[int] = None,
250 x_min: Optional[float] = None,
251 x_max: Optional[float] = None,
252 y_min: Optional[float] = None,
253 y_max: Optional[float] = None,
254 top: Optional[GridDataArray] = None,
255 bottom: Optional[GridDataArray] = None,
256 state_for_boundary: Optional[GridDataArray] = None,
257 ) -> Package:
258 """
259 Clip a package by a bounding box (time, layer, y, x).
261 The well package doesn't use the layer attribute to describe its depth and length.
262 Instead, it uses the screen_top and screen_bottom parameters which corresponds with
263 the z-coordinates of the top and bottom of the well. To go from a layer_min and
264 layer_max to z-values used for clipping the well a top and bottom array have to be
265 provided as well.
267 Slicing intervals may be half-bounded, by providing None:
269 * To select 500.0 <= x <= 1000.0:
270 ``clip_box(x_min=500.0, x_max=1000.0)``.
271 * To select x <= 1000.0: ``clip_box(x_min=None, x_max=1000.0)``
272 or ``clip_box(x_max=1000.0)``.
273 * To select x >= 500.0: ``clip_box(x_min = 500.0, x_max=None.0)``
274 or ``clip_box(x_min=1000.0)``.
276 Parameters
277 ----------
278 time_min: optional
279 time_max: optional
280 layer_min: optional, int
281 layer_max: optional, int
282 x_min: optional, float
283 x_max: optional, float
284 y_min: optional, float
285 y_max: optional, float
286 top: optional, GridDataArray
287 bottom: optional, GridDataArray
288 state_for_boundary: optional, GridDataArray
290 Returns
291 -------
292 sliced : Package
293 """
294 if (layer_max or layer_min) and (top is None or bottom is None):
295 raise ValueError(
296 "When clipping by layer both the top and bottom should be defined"
297 )
299 if top is not None:
300 # Bug in mypy when using unions in isInstance
301 if not isinstance(top, GridDataArray) or "layer" not in top.coords: # type: ignore
302 top = create_layered_top(bottom, top)
304 # The super method will select in the time dimension without issues.
305 new = super().clip_box(time_min=time_min, time_max=time_max)
307 ds = new.dataset
309 z_max = self._find_well_value_at_layer(ds, top, layer_max)
310 z_min = self._find_well_value_at_layer(ds, bottom, layer_min)
312 if z_max is not None:
313 ds["screen_top"] = ds["screen_top"].clip(None, z_max)
314 if z_min is not None:
315 ds["screen_bottom"] = ds["screen_bottom"].clip(z_min, None)
317 # Initiate array of True with right shape to deal with case no spatial
318 # selection needs to be done.
319 in_bounds = np.full(ds.dims["index"], True)
320 # Select all variables along "index" dimension
321 in_bounds &= values_within_range(ds["x"], x_min, x_max)
322 in_bounds &= values_within_range(ds["y"], y_min, y_max)
323 in_bounds &= values_within_range(ds["screen_top"], z_min, z_max)
324 in_bounds &= values_within_range(ds["screen_bottom"], z_min, z_max)
325 # remove wells where the screen bottom and top are the same
326 in_bounds &= abs(ds["screen_bottom"] - ds["screen_top"]) > 1e-5
327 # Replace dataset with reduced dataset based on booleans
328 new.dataset = ds.loc[{"index": in_bounds}]
330 return new
332 @staticmethod
333 def _find_well_value_at_layer(
334 well_dataset: xr.Dataset, grid: GridDataArray, layer: Optional[int]
335 ):
336 value = None if layer is None else grid.isel(layer=layer)
338 # if value is a grid select the values at the well locations and drop the dimensions
339 if (value is not None) and is_spatial_2D(value):
340 value = imod.select.points_values(
341 value,
342 x=well_dataset["x"].values,
343 y=well_dataset["y"].values,
344 out_of_bounds="ignore",
345 ).drop_vars(lambda x: x.coords)
347 return value
349 def write(
350 self,
351 pkgname: str,
352 globaltimes: Union[list[np.datetime64], np.ndarray],
353 write_context: WriteContext,
354 ):
355 raise NotImplementedError(
356 "To write a wel package first convert it to a MF6 well using to_mf6_pkg."
357 )
359 def __create_wells_df(self) -> pd.DataFrame:
360 wells_df = self.dataset.to_dataframe()
361 wells_df = wells_df.rename(
362 columns={
363 "screen_top": "top",
364 "screen_bottom": "bottom",
365 }
366 )
368 return wells_df
370 def __create_assigned_wells(
371 self,
372 wells_df: pd.DataFrame,
373 active: GridDataArray,
374 top: GridDataArray,
375 bottom: GridDataArray,
376 k: GridDataArray,
377 minimum_k: float,
378 minimum_thickness: float,
379 ):
380 # Ensure top, bottom & k
381 # are broadcasted to 3d grid
382 like = ones_like(active)
383 bottom = like * bottom
384 top_2d = (like * top).sel(layer=1)
385 top_3d = bottom.shift(layer=1).fillna(top_2d)
387 k = like * k
389 index_names = wells_df.index.names
391 # Unset multi-index, because assign_wells cannot deal with
392 # multi-indices which is returned by self.dataset.to_dataframe() in
393 # case of a "time" and "species" coordinate.
394 wells_df = wells_df.reset_index()
396 wells_assigned = assign_wells(
397 wells_df, top_3d, bottom, k, minimum_thickness, minimum_k
398 )
399 # Set multi-index again
400 wells_assigned = wells_assigned.set_index(index_names).sort_index()
402 return wells_assigned
404 def __create_dataset_vars(
405 self, wells_assigned: pd.DataFrame, wells_df: pd.DataFrame, cellid: xr.DataArray
406 ) -> xr.Dataset:
407 """
408 Create dataset with all variables (rate, concentration), with a similar shape as the cellids.
409 """
410 data_vars = ["rate"]
411 if "concentration" in wells_assigned.columns:
412 data_vars.append("concentration")
414 ds_vars = wells_assigned[data_vars].to_xarray()
415 # "rate" variable in conversion from multi-indexed DataFrame to xarray
416 # DataArray results in duplicated values for "rate" along dimension
417 # "species". Select first species to reduce this again.
418 index_names = wells_df.index.names
419 if "species" in index_names:
420 ds_vars["rate"] = ds_vars["rate"].isel(species=0)
422 # Carefully rename the dimension and set coordinates
423 d_rename = {"index": "ncellid"}
424 ds_vars = ds_vars.rename_dims(**d_rename).rename_vars(**d_rename)
425 ds_vars = ds_vars.assign_coords(**{"ncellid": cellid.coords["ncellid"].values})
427 return ds_vars
429 def __create_cellid(self, wells_assigned: pd.DataFrame, active: xr.DataArray):
430 like = ones_like(active)
432 # Groupby index and select first, to unset any duplicate records
433 # introduced by the multi-indexed "time" dimension.
434 df_for_cellid = wells_assigned.groupby("index").first()
435 d_for_cellid = df_for_cellid[["x", "y", "layer"]].to_dict("list")
437 return self.__derive_cellid_from_points(like, **d_for_cellid)
439 @staticmethod
440 def __derive_cellid_from_points(
441 dst_grid: GridDataArray,
442 x: list,
443 y: list,
444 layer: list,
445 ) -> GridDataArray:
446 """
447 Create DataArray with Modflow6 cell identifiers based on x, y coordinates
448 in a dataframe. For structured grid this DataArray contains 3 columns:
449 ``layer, row, column``. For unstructured grids, this contains 2 columns:
450 ``layer, cell2d``.
451 See also: https://water.usgs.gov/water-resources/software/MODFLOW-6/mf6io_6.4.0.pdf#page=35
453 Note
454 ----
455 The "layer" coordinate should already be provided in the dataframe.
456 To determine the layer coordinate based on screen depts, look at
457 :func:`imod.prepare.wells.assign_wells`.
459 Parameters
460 ----------
461 dst_grid: {xr.DataArray, xu.UgridDataArray}
462 Destination grid to map the points to based on their x and y coordinates.
463 x: {list, np.array}
464 array-like with x-coordinates
465 y: {list, np.array}
466 array-like with y-coordinates
467 layer: {list, np.array}
468 array-like with layer-coordinates
470 Returns
471 -------
472 cellid : xr.DataArray
473 2D DataArray with a ``ncellid`` rows and 3 to 2 columns, depending
474 on whether on a structured or unstructured grid."""
476 # Find indices belonging to x, y coordinates
477 indices_cell2d = points_indices(dst_grid, out_of_bounds="ignore", x=x, y=y)
478 # Convert cell2d indices from 0-based to 1-based.
479 indices_cell2d = dict((dim, index + 1) for dim, index in indices_cell2d.items())
480 # Prepare layer indices, for later concatenation
482 if isinstance(dst_grid, xu.UgridDataArray):
483 indices_layer = xr.DataArray(
484 layer, coords=indices_cell2d["mesh2d_nFaces"].coords
485 )
486 face_dim = dst_grid.ugrid.grid.face_dimension
487 indices_cell2d_dims = [face_dim]
488 cell2d_coords = ["cell2d"]
489 else:
490 indices_layer = xr.DataArray(layer, coords=indices_cell2d["x"].coords)
491 indices_cell2d_dims = ["y", "x"]
492 cell2d_coords = ["row", "column"]
494 # Prepare cellid array of the right shape.
495 cellid_ls = [indices_layer] + [
496 indices_cell2d[dim] for dim in indices_cell2d_dims
497 ]
498 cellid = xr.concat(cellid_ls, dim="nmax_cellid")
499 # Rename generic dimension name "index" to ncellid.
500 cellid = cellid.rename(index="ncellid")
501 # Put dimensions in right order after concatenation.
502 cellid = cellid.transpose("ncellid", "nmax_cellid")
503 # Assign extra coordinate names.
504 coords = {
505 "nmax_cellid": ["layer"] + cell2d_coords,
506 "x": ("ncellid", x),
507 "y": ("ncellid", y),
508 }
509 cellid = cellid.assign_coords(**coords)
511 return cellid
513 def render(self, directory, pkgname, globaltimes, binary):
514 raise NotImplementedError(
515 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()"
516 )
518 def to_mf6_pkg(
519 self,
520 active: GridDataArray,
521 top: GridDataArray,
522 bottom: GridDataArray,
523 k: GridDataArray,
524 validate: bool = False,
525 is_partitioned: bool = False,
526 ) -> Mf6Wel:
527 """
528 Write package to Modflow 6 package.
530 Based on the model grid and top and bottoms, cellids are determined.
531 When well screens hit multiple layers, groundwater extractions are
532 distributed based on layer transmissivities. Wells located in inactive
533 cells are removed.
535 Note
536 ----
537 The well distribution based on transmissivities assumes confined
538 aquifers. If wells fall dry (and the rate distribution has to be
539 recomputed at runtime), it is better to use the Multi-Aquifer Well
540 package.
542 Parameters
543 ----------
544 is_partitioned: bool
545 validate: bool
546 Run validation before converting
547 active: {xarry.DataArray, xugrid.UgridDataArray}
548 Grid with active cells.
549 top: {xarry.DataArray, xugrid.UgridDataArray}
550 Grid with top of model layers.
551 bottom: {xarry.DataArray, xugrid.UgridDataArray}
552 Grid with bottom of model layers.
553 k: {xarry.DataArray, xugrid.UgridDataArray}
554 Grid with hydraulic conductivities.
555 Returns
556 -------
557 Mf6Wel
558 Object with wells as list based input.
559 """
560 if validate:
561 errors = self._validate(self._write_schemata)
562 if len(errors) > 0:
563 message = validation_pkg_error_message(errors)
564 raise ValidationError(message)
566 minimum_k = self.dataset["minimum_k"].item()
567 minimum_thickness = self.dataset["minimum_thickness"].item()
569 wells_df = self.__create_wells_df()
570 wells_assigned = self.__create_assigned_wells(
571 wells_df, active, top, bottom, k, minimum_k, minimum_thickness
572 )
574 nwells_df = len(wells_df["id"].unique())
575 nwells_assigned = (
576 0 if wells_assigned.empty else len(wells_assigned["id"].unique())
577 )
579 if nwells_df == 0:
580 raise ValueError("No wells were assigned in package. None were present.")
581 # @TODO: reinstate this check. issue github #621.
582 if not is_partitioned and nwells_df != nwells_assigned:
583 raise ValueError(
584 "One or more well(s) are completely invalid due to minimum conductivity and thickness constraints."
585 )
587 ds = xr.Dataset()
588 ds["cellid"] = self.__create_cellid(wells_assigned, active)
590 ds_vars = self.__create_dataset_vars(wells_assigned, wells_df, ds["cellid"])
591 ds = ds.assign(**ds_vars.data_vars)
593 ds = remove_inactive(ds, active)
594 ds["save_flows"] = self["save_flows"].values[()]
595 ds["print_flows"] = self["print_flows"].values[()]
596 ds["print_input"] = self["print_input"].values[()]
598 return Mf6Wel(**ds.data_vars)
600 def mask(self, domain: GridDataArray) -> Well:
601 """
602 Mask wells based on two-dimensional domain. For three-dimensional
603 masking: Wells falling in inactive cells are automatically removed in
604 the call to write to Modflow 6 package. You can verify this by calling
605 the ``to_mf6_pkg`` method.
606 """
608 # Drop layer coordinate if present, otherwise a layer coordinate is assigned
609 # which causes conflicts downstream when assigning wells and deriving
610 # cellids.
611 domain_2d = domain.isel(layer=0, drop=True, missing_dims="ignore").drop_vars(
612 "layer", errors="ignore"
613 )
614 return mask_2D(self, domain_2d)
616 def get_regrid_methods(self) -> Optional[dict[str, Tuple[RegridderType, str]]]:
617 return self._regrid_method
620class WellDisStructured(DisStructuredBoundaryCondition):
621 """
622 WEL package for structured discretization (DIS) models .
623 Any number of WEL Packages can be specified for a single groundwater flow model.
624 https://water.usgs.gov/water-resources/software/MODFLOW-6/mf6io_6.0.4.pdf#page=63
626 .. warning::
627 This class is deprecated and will be deleted in a future release.
628 Consider changing your code to use the ``imod.mf6.Well`` package.
630 Parameters
631 ----------
632 layer: list of int
633 Model layer in which the well is located.
634 row: list of int
635 Row in which the well is located.
636 column: list of int
637 Column in which the well is located.
638 rate: float or list of floats
639 is the volumetric well rate. A positive value indicates well
640 (injection) and a negative value indicates discharge (extraction) (q).
641 concentration: array of floats (xr.DataArray, optional)
642 if this flow package is used in simulations also involving transport, then this array is used
643 as the concentration for inflow over this boundary.
644 concentration_boundary_type: ({"AUX", "AUXMIXED"}, optional)
645 if this flow package is used in simulations also involving transport, then this keyword specifies
646 how outflow over this boundary is computed.
647 print_input: ({True, False}, optional)
648 keyword to indicate that the list of well information will be written to
649 the listing file immediately after it is read.
650 Default is False.
651 print_flows: ({True, False}, optional)
652 Indicates that the list of well flow rates will be printed to the
653 listing file for every stress period time step in which "BUDGET PRINT"
654 is specified in Output Control. If there is no Output Control option
655 and PRINT FLOWS is specified, then flow rates are printed for the last
656 time step of each stress period.
657 Default is False.
658 save_flows: ({True, False}, optional)
659 Indicates that well flow terms will be written to the file specified
660 with "BUDGET FILEOUT" in Output Control.
661 Default is False.
662 observations: [Not yet supported.]
663 Default is None.
664 validate: {True, False}
665 Flag to indicate whether the package should be validated upon
666 initialization. This raises a ValidationError if package input is
667 provided in the wrong manner. Defaults to True.
668 repeat_stress: Optional[xr.DataArray] of datetimes
669 Used to repeat data for e.g. repeating stress periods such as
670 seasonality without duplicating the values. The DataArray should have
671 dimensions ``("repeat", "repeat_items")``. The ``repeat_items``
672 dimension should have size 2: the first value is the "key", the second
673 value is the "value". For the "key" datetime, the data of the "value"
674 datetime will be used. Can also be set with a dictionary using the
675 ``set_repeat_stress`` method.
676 """
678 _pkg_id = "wel"
679 _period_data = ("layer", "row", "column", "rate")
680 _keyword_map = {}
681 _template = DisStructuredBoundaryCondition._initialize_template(_pkg_id)
682 _auxiliary_data = {"concentration": "species"}
684 _init_schemata = {
685 "layer": [DTypeSchema(np.integer)],
686 "row": [DTypeSchema(np.integer)],
687 "column": [DTypeSchema(np.integer)],
688 "rate": [DTypeSchema(np.floating)],
689 "concentration": [DTypeSchema(np.floating)],
690 }
692 _write_schemata = {}
694 @init_log_decorator()
695 def __init__(
696 self,
697 layer,
698 row,
699 column,
700 rate,
701 concentration=None,
702 concentration_boundary_type="aux",
703 print_input=False,
704 print_flows=False,
705 save_flows=False,
706 observations=None,
707 validate: bool = True,
708 repeat_stress=None,
709 ):
710 dict_dataset = {
711 "layer": _assign_dims(layer),
712 "row": _assign_dims(row),
713 "column": _assign_dims(column),
714 "rate": _assign_dims(rate),
715 "print_input": print_input,
716 "print_flows": print_flows,
717 "save_flows": save_flows,
718 "observations": observations,
719 "repeat_stress": repeat_stress,
720 "concentration": concentration,
721 "concentration_boundary_type": concentration_boundary_type,
722 }
723 super().__init__(dict_dataset)
724 self._validate_init_schemata(validate)
726 warnings.warn(
727 f"{self.__class__.__name__} is deprecated and will be removed in the v1.0 release."
728 "Please adapt your code to use the imod.mf6.Well package",
729 DeprecationWarning,
730 )
732 def clip_box(
733 self,
734 time_min: Optional[cftime.datetime | np.datetime64 | str] = None,
735 time_max: Optional[cftime.datetime | np.datetime64 | str] = None,
736 layer_min: Optional[int] = None,
737 layer_max: Optional[int] = None,
738 x_min: Optional[float] = None,
739 x_max: Optional[float] = None,
740 y_min: Optional[float] = None,
741 y_max: Optional[float] = None,
742 top: Optional[GridDataArray] = None,
743 bottom: Optional[GridDataArray] = None,
744 state_for_boundary: Optional[GridDataArray] = None,
745 ) -> Package:
746 """
747 Clip a package by a bounding box (time, layer, y, x).
749 Slicing intervals may be half-bounded, by providing None:
751 * To select 500.0 <= x <= 1000.0:
752 ``clip_box(x_min=500.0, x_max=1000.0)``.
753 * To select x <= 1000.0: ``clip_box(x_min=None, x_max=1000.0)``
754 or ``clip_box(x_max=1000.0)``.
755 * To select x >= 500.0: ``clip_box(x_min = 500.0, x_max=None.0)``
756 or ``clip_box(x_min=1000.0)``.
758 Parameters
759 ----------
760 time_min: optional
761 time_max: optional
762 layer_min: optional, int
763 layer_max: optional, int
764 x_min: optional, float
765 x_max: optional, float
766 y_min: optional, float
767 y_max: optional, float
768 top: optional, GridDataArray
769 bottom: optional, GridDataArray
770 state_for_boundary: optional, GridDataArray
772 Returns
773 -------
774 sliced : Package
775 """
776 # TODO: include x and y values.
777 for arg in (
778 layer_min,
779 layer_max,
780 x_min,
781 x_max,
782 y_min,
783 y_max,
784 ):
785 if arg is not None:
786 raise NotImplementedError("Can only clip_box in time for Well packages")
788 # The super method will select in the time dimension without issues.
789 new = super().clip_box(time_min=time_min, time_max=time_max)
790 return new
793class WellDisVertices(DisVerticesBoundaryCondition):
794 """
795 WEL package for discretization by vertices (DISV) models. Any number of WEL
796 Packages can be specified for a single groundwater flow model.
797 https://water.usgs.gov/water-resources/software/MODFLOW-6/mf6io_6.0.4.pdf#page=63
799 .. warning::
800 This class is deprecated and will be deleted in a future release.
801 Consider changing your code to use the ``imod.mf6.Well`` package.
803 Parameters
804 ----------
805 layer: list of int
806 Modellayer in which the well is located.
807 cell2d: list of int
808 Cell in which the well is located.
809 rate: float or list of floats
810 is the volumetric well rate. A positive value indicates well (injection)
811 and a negative value indicates discharge (extraction) (q).
812 concentration: array of floats (xr.DataArray, optional)
813 if this flow package is used in simulations also involving transport,
814 then this array is used as the concentration for inflow over this
815 boundary.
816 concentration_boundary_type: ({"AUX", "AUXMIXED"}, optional)
817 if this flow package is used in simulations also involving transport,
818 then this keyword specifies how outflow over this boundary is computed.
819 print_input: ({True, False}, optional)
820 keyword to indicate that the list of well information will be written to
821 the listing file immediately after it is read. Default is False.
822 print_flows: ({True, False}, optional)
823 Indicates that the list of well flow rates will be printed to the
824 listing file for every stress period time step in which "BUDGET PRINT"
825 is specified in Output Control. If there is no Output Control option and
826 PRINT FLOWS is specified, then flow rates are printed for the last time
827 step of each stress period. Default is False.
828 save_flows: ({True, False}, optional)
829 Indicates that well flow terms will be written to the file specified
830 with "BUDGET FILEOUT" in Output Control. Default is False.
831 observations: [Not yet supported.]
832 Default is None.
833 validate: {True, False}
834 Flag to indicate whether the package should be validated upon
835 initialization. This raises a ValidationError if package input is
836 provided in the wrong manner. Defaults to True.
837 """
839 _pkg_id = "wel"
840 _period_data = ("layer", "cell2d", "rate")
841 _keyword_map = {}
842 _template = DisVerticesBoundaryCondition._initialize_template(_pkg_id)
843 _auxiliary_data = {"concentration": "species"}
845 _init_schemata = {
846 "layer": [DTypeSchema(np.integer)],
847 "cell2d": [DTypeSchema(np.integer)],
848 "rate": [DTypeSchema(np.floating)],
849 "concentration": [DTypeSchema(np.floating)],
850 }
852 _write_schemata = {}
854 @init_log_decorator()
855 def __init__(
856 self,
857 layer,
858 cell2d,
859 rate,
860 concentration=None,
861 concentration_boundary_type="aux",
862 print_input=False,
863 print_flows=False,
864 save_flows=False,
865 observations=None,
866 validate: bool = True,
867 ):
868 dict_dataset = {
869 "layer": _assign_dims(layer),
870 "cell2d": _assign_dims(cell2d),
871 "rate": _assign_dims(rate),
872 "print_input": print_input,
873 "print_flows": print_flows,
874 "save_flows": save_flows,
875 "observations": observations,
876 "concentration": concentration,
877 "concentration_boundary_type": concentration_boundary_type,
878 }
879 super().__init__(dict_dataset)
880 self._validate_init_schemata(validate)
882 warnings.warn(
883 f"{self.__class__.__name__} is deprecated and will be removed in the v1.0 release."
884 "Please adapt your code to use the imod.mf6.Well package",
885 DeprecationWarning,
886 )
888 def clip_box(
889 self,
890 time_min: Optional[cftime.datetime | np.datetime64 | str] = None,
891 time_max: Optional[cftime.datetime | np.datetime64 | str] = None,
892 layer_min: Optional[int] = None,
893 layer_max: Optional[int] = None,
894 x_min: Optional[float] = None,
895 x_max: Optional[float] = None,
896 y_min: Optional[float] = None,
897 y_max: Optional[float] = None,
898 top: Optional[GridDataArray] = None,
899 bottom: Optional[GridDataArray] = None,
900 state_for_boundary: Optional[GridDataArray] = None,
901 ) -> Package:
902 """
903 Clip a package by a bounding box (time, layer, y, x).
905 Slicing intervals may be half-bounded, by providing None:
907 * To select 500.0 <= x <= 1000.0:
908 ``clip_box(x_min=500.0, x_max=1000.0)``.
909 * To select x <= 1000.0: ``clip_box(x_min=None, x_max=1000.0)``
910 or ``clip_box(x_max=1000.0)``.
911 * To select x >= 500.0: ``clip_box(x_min = 500.0, x_max=None.0)``
912 or ``clip_box(x_min=1000.0)``.
914 Parameters
915 ----------
916 time_min: optional
917 time_max: optional
918 layer_min: optional, int
919 layer_max: optional, int
920 x_min: optional, float
921 x_max: optional, float
922 y_min: optional, float
923 y_max: optional, float
924 top: optional, GridDataArray
925 bottom: optional, GridDataArray
926 state_for_boundary: optional, GridDataArray
928 Returns
929 -------
930 clipped: Package
931 """
932 # TODO: include x and y values.
933 for arg in (
934 layer_min,
935 layer_max,
936 x_min,
937 x_max,
938 y_min,
939 y_max,
940 ):
941 if arg is not None:
942 raise NotImplementedError("Can only clip_box in time for Well packages")
944 # The super method will select in the time dimension without issues.
945 new = super().clip_box(time_min=time_min, time_max=time_max)
946 return new