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

1from __future__ import annotations 

2 

3import warnings 

4from typing import Any, Optional, Tuple, Union 

5 

6import cftime 

7import numpy as np 

8import numpy.typing as npt 

9import pandas as pd 

10import xarray as xr 

11import xugrid as xu 

12 

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 

39 

40 

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 

56 

57 

58def mask_2D(package: Well, domain_2d: GridDataArray) -> Well: 

59 point_active = points_values(domain_2d, x=package.x, y=package.y) 

60 

61 is_inside_exterior = point_active == 1 

62 selection = package.dataset.loc[{"index": is_inside_exterior}] 

63 

64 cls = type(package) 

65 new = cls.__new__(cls) 

66 new.dataset = selection 

67 return new 

68 

69 

70class Well(BoundaryCondition, IPointDataPackage): 

71 """ 

72 Agnostic WEL package, which accepts x, y and a top and bottom of the well screens. 

73 

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 

77 

78 Parameters 

79 ---------- 

80 

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. 

138 

139 Examples 

140 --------- 

141 

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] 

147 

148 >>> imod.mf6.Well(x, y, screen_top, screen_bottom, rate) 

149 

150 For a transient well: 

151 

152 >>> weltimes = pd.date_range("2000-01-01", "2000-01-03") 

153 

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",)) 

156 

157 >>> imod.mf6.Well(x, y, screen_top, screen_bottom, rate_transient) 

158 """ 

159 

160 @property 

161 def x(self) -> npt.NDArray[np.float64]: 

162 return self.dataset["x"].values 

163 

164 @property 

165 def y(self) -> npt.NDArray[np.float64]: 

166 return self.dataset["y"].values 

167 

168 _pkg_id = "wel" 

169 

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 } 

187 

188 _regrid_method: dict[str, Tuple[RegridderType, str]] = {} 

189 

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) 

239 

240 @classmethod 

241 def is_grid_agnostic_package(cls) -> bool: 

242 return True 

243 

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). 

260 

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. 

266 

267 Slicing intervals may be half-bounded, by providing None: 

268 

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)``. 

275 

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 

289 

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 ) 

298 

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) 

303 

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) 

306 

307 ds = new.dataset 

308 

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) 

311 

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) 

316 

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}] 

329 

330 return new 

331 

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) 

337 

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) 

346 

347 return value 

348 

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 ) 

358 

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 ) 

367 

368 return wells_df 

369 

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) 

386 

387 k = like * k 

388 

389 index_names = wells_df.index.names 

390 

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() 

395 

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() 

401 

402 return wells_assigned 

403 

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") 

413 

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) 

421 

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}) 

426 

427 return ds_vars 

428 

429 def __create_cellid(self, wells_assigned: pd.DataFrame, active: xr.DataArray): 

430 like = ones_like(active) 

431 

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") 

436 

437 return self.__derive_cellid_from_points(like, **d_for_cellid) 

438 

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 

452 

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`. 

458 

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 

469 

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.""" 

475 

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 

481 

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"] 

493 

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) 

510 

511 return cellid 

512 

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 ) 

517 

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. 

529 

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. 

534 

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. 

541 

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) 

565 

566 minimum_k = self.dataset["minimum_k"].item() 

567 minimum_thickness = self.dataset["minimum_thickness"].item() 

568 

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 ) 

573 

574 nwells_df = len(wells_df["id"].unique()) 

575 nwells_assigned = ( 

576 0 if wells_assigned.empty else len(wells_assigned["id"].unique()) 

577 ) 

578 

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 ) 

586 

587 ds = xr.Dataset() 

588 ds["cellid"] = self.__create_cellid(wells_assigned, active) 

589 

590 ds_vars = self.__create_dataset_vars(wells_assigned, wells_df, ds["cellid"]) 

591 ds = ds.assign(**ds_vars.data_vars) 

592 

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[()] 

597 

598 return Mf6Wel(**ds.data_vars) 

599 

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 """ 

607 

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) 

615 

616 def get_regrid_methods(self) -> Optional[dict[str, Tuple[RegridderType, str]]]: 

617 return self._regrid_method 

618 

619 

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 

625 

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. 

629 

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 """ 

677 

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"} 

683 

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 } 

691 

692 _write_schemata = {} 

693 

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) 

725 

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 ) 

731 

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). 

748 

749 Slicing intervals may be half-bounded, by providing None: 

750 

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)``. 

757 

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 

771 

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") 

787 

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 

791 

792 

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 

798 

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. 

802 

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 """ 

838 

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"} 

844 

845 _init_schemata = { 

846 "layer": [DTypeSchema(np.integer)], 

847 "cell2d": [DTypeSchema(np.integer)], 

848 "rate": [DTypeSchema(np.floating)], 

849 "concentration": [DTypeSchema(np.floating)], 

850 } 

851 

852 _write_schemata = {} 

853 

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) 

881 

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 ) 

887 

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). 

904 

905 Slicing intervals may be half-bounded, by providing None: 

906 

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)``. 

913 

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 

927 

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") 

943 

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