Coverage for C:\src\imod-python\imod\mf6\npf.py: 89%

53 statements  

« prev     ^ index     » next       coverage.py v7.5.1, created at 2024-05-08 14:15 +0200

1import warnings 

2from typing import Optional, Tuple 

3 

4import numpy as np 

5 

6from imod.logging import init_log_decorator 

7from imod.mf6.interfaces.iregridpackage import IRegridPackage 

8from imod.mf6.package import Package 

9from imod.mf6.utilities.regrid import RegridderType 

10from imod.mf6.validation import PKG_DIMS_SCHEMA 

11from imod.schemata import ( 

12 AllValueSchema, 

13 CompatibleSettingsSchema, 

14 DimsSchema, 

15 DTypeSchema, 

16 IdentityNoDataSchema, 

17 IndexesSchema, 

18) 

19from imod.typing import GridDataArray 

20 

21 

22def _dataarray_to_bool(griddataarray: GridDataArray) -> bool: 

23 if griddataarray is None or griddataarray.values is None: 

24 return False 

25 

26 if griddataarray.values.size != 1: 

27 raise ValueError("DataArray is not a single value") 

28 

29 if griddataarray.values.dtype != bool: 

30 raise ValueError("DataArray is not a boolean") 

31 

32 return griddataarray.values.item() 

33 

34 

35class NodePropertyFlow(Package, IRegridPackage): 

36 """ 

37 Node Property Flow package. 

38 

39 In this package the hydraulic conductivity and rewetting in the model is 

40 specified. A single NPF Package is required for each GWF model. 

41 https://water.usgs.gov/water-resources/software/MODFLOW-6/mf6io_6.0.4.pdf#page=51 

42 

43 A note about regridding: the fields k, k22, k33 define the principal 

44 components of an anisotropic conductivity tensor. By default, k and k22 are 

45 in the horizontal plane and k33 is vertical. Angle1, angle2 and angle3 

46 define the rotation of this tensor. The regridding methods associated by 

47 default are chosen based on the assumption that k and k22 are horizontal and 

48 k33 is vertical. If this is not the case, it is up to the user to regrid the 

49 npf package using other regridding methods. This may be recommended if for 

50 example the rotation is such that k has become vertical and k33 horizontal. 

51 

52 Parameters 

53 ---------- 

54 icelltype: array of int (xr.DataArray) 

55 flag for each cell that specifies how saturated thickness is treated. 0 

56 means saturated thickness is held constant; >0 means saturated thickness 

57 varies with computed head when head is below the cell top; <0 means 

58 saturated thickness varies with computed head unless the 

59 starting_head_as_confined_thickness option is in effect. When 

60 starting_head_as_confined_thickness is in effect, a negative value of 

61 icelltype indicates that saturated thickness will be computed as 

62 strt-bot and held constant. 

63 k: array of floats (xr.DataArray) 

64 is the hydraulic conductivity. For the common case in which the user 

65 would like to specify the horizontal hydraulic conductivity and the 

66 vertical hydraulic conductivity, then K should be assigned as the 

67 horizontal hydraulic conductivity, K33 should be assigned as the 

68 vertical hydraulic conductivity, and K22 and the three rotation 

69 angles should not be specified. When more sophisticated anisotropy is 

70 required, then K corresponds to the K11 hydraulic conductivity axis. All 

71 included cells (idomain > 0) must have a K value greater than zero 

72 rewet: ({True, False}, optional) 

73 activates model rewetting. 

74 Default is False. 

75 rewet_layer: float 

76 is a combination of the wetting threshold and a flag to indicate which 

77 neighboring cells can cause a cell to become wet. If rewet_layer < 0, 

78 only a cell below a dry cell can cause the cell to become wet. If 

79 rewet_layer > 0, the cell below a dry cell and horizontally adjacent 

80 cells can cause a cell to become wet. If rewet_layer is 0, the cell 

81 cannot be wetted. The absolute value of rewet_layer is the wetting 

82 threshold. When the sum of BOT and the absolute value of rewet_layer at 

83 a dry cell is equaled or exceeded by the head at an adjacent cell, the 

84 cell is wetted. rewet_layer must be specified if "rewet" is specified in 

85 the OPTIONS block. If "rewet" is not specified in the options block, 

86 then rewet_layer can be entered, and memory will be allocated for it, 

87 even though it is not used. (WETDRY) 

88 Default is None. 

89 rewet_factor: 

90 is a keyword and factor that is included in the calculation of the head 

91 that is initially established at a cell when that cell is converted from 

92 dry to wet. (WETFCT) 

93 Default is None. 

94 rewet_iterations: 

95 is a keyword and iteration interval for attempting to wet cells. Wetting 

96 is attempted every rewet_iterations iteration. This applies to outer 

97 iterations and not inner iterations. If rewet_iterations is specified as 

98 zero or less, then the value is changed to 1. (IWETIT) 

99 Default is None. 

100 rewet_method: 

101 is a keyword and integer flag that determines which equation is used to 

102 define the initial head at cells that become wet. If rewet_method is 0, 

103 h = BOT + rewet_factor (hm - BOT). If rewet_method is not 0, h = BOT + 

104 rewet_factor (THRESH). (IHDWET) 

105 Default is None. 

106 k22: array of floats (xr.DataArray) 

107 is the hydraulic conductivity of the second ellipsoid axis; for an 

108 unrotated case this is the hydraulic conductivity in the y direction. If 

109 K22 is not included, then K22 is set equal to K. 

110 For a regular MODFLOW grid (DIS Package is used) in which no rotation 

111 angles are specified, K22 is the hydraulic conductivity along columns in 

112 the y direction. For an unstructured DISU grid, the user must assign 

113 principal x and y axes and provide the angle for each cell face relative 

114 to the assigned x direction. All included cells (idomain > 0) must have 

115 a K22 value greater than zero. 

116 Default is None. 

117 k33: array of floats (xr.DataArray) 

118 is the hydraulic conductivity of the third ellipsoid axis; for an 

119 unrotated case, this is the vertical hydraulic conductivity. When 

120 anisotropy is applied, K33 corresponds to the K33 tensor component. All 

121 included cells (idomain > 0) must have a K33 value greater than zero. 

122 Default is None. 

123 angle1: float 

124 is a rotation angle of the hydraulic conductivity tensor in degrees. The 

125 angle represents the first of three sequential rotations of the 

126 hydraulic conductivity ellipsoid. With the K11, K22, and K33 axes of the 

127 ellipsoid initially aligned with the x, y, and z coordinate axes, 

128 respectively, angle1 rotates the ellipsoid about its K33 axis (within 

129 the x - y plane). A positive value represents counter-clockwise rotation 

130 when viewed from any point on the positive K33 axis, looking toward the 

131 center of the ellipsoid. A value of zero indicates that the K11 axis 

132 lies within the x - z plane. If angle1 is not specified, default values 

133 of zero are assigned to angle1, angle2, and angle3, in which case the 

134 K11, K22, and K33 axes are aligned with the x, y, and z axes, 

135 respectively. 

136 Default is None. 

137 angle2: float 

138 is a rotation angle of the hydraulic conductivity tensor in degrees. The 

139 angle represents the second of three sequential rotations of the 

140 hydraulic conductivity ellipsoid. Following the rotation by angle1 

141 described above, angle2 rotates the ellipsoid about its K22 axis (out of 

142 the x - y plane). An array can be specified for angle2 only if angle1 is 

143 also specified. A positive value of angle2 represents clockwise rotation 

144 when viewed from any point on the positive K22 axis, looking toward the 

145 center of the ellipsoid. A value of zero indicates that the K11 axis 

146 lies within the x - y plane. If angle2 is not specified, default values 

147 of zero are assigned to angle2 and angle3; connections that are not 

148 user-designated as vertical are assumed to be strictly horizontal (that 

149 is, to have no z component to their orientation); and connection lengths 

150 are based on horizontal distances. 

151 Default is None. 

152 angle3: float 

153 is a rotation angle of the hydraulic conductivity tensor in degrees. The 

154 angle represents the third of three sequential rotations of the 

155 hydraulic conductivity ellipsoid. Following the rotations by angle1 and 

156 angle2 described above, angle3 rotates the ellipsoid about its K11 axis. 

157 An array can be specified for angle3 only if angle1 and angle2 are also 

158 specified. An array must be specified for angle3 if angle2 is specified. 

159 A positive value of angle3 represents clockwise rotation when viewed 

160 from any point on the positive K11 axis, looking toward the center of 

161 the ellipsoid. A value of zero indicates that the K22 axis lies within 

162 the x - y plane. 

163 Default is None. 

164 alternative_cell_averaging : str 

165 Method calculating horizontal cell connection conductance. 

166 Options: {"LOGARITHMIC", "AMT-LMK", or "AMT-HMK"} 

167 Default: uses harmonic mean for averaging 

168 save_flows: ({True, False}, optional) 

169 keyword to indicate that cell-by-cell flow terms will be written to the 

170 file specified with "budget save file" in Output Control. 

171 Default is False. 

172 starting_head_as_confined_thickness: ({True, False}, optional) 

173 indicates that cells having a negative icelltype are confined, and their 

174 cell thickness for conductance calculations will be computed as strt-bot 

175 rather than top-bot. 

176 (THICKSTRT) 

177 Default is False. 

178 variable_vertical_conductance: ({True, False}, optional) 

179 keyword to indicate that the vertical conductance will be calculated 

180 using the saturated thickness and properties of the overlying cell and 

181 the thickness and properties of the underlying cell. if the dewatered 

182 keyword is also specified, then the vertical conductance is calculated 

183 using only the saturated thickness and properties of the overlying cell 

184 if the head in the underlying cell is below its top. if these keywords 

185 are not specified, then the default condition is to calculate the 

186 vertical conductance at the start of the simulation using the initial 

187 head and the cell properties. the vertical conductance remains constant 

188 for the entire simulation. 

189 (VARIABLECV) 

190 Default is False. 

191 dewatered: ({True, False}, optional) 

192 If the dewatered keyword is specified, then the vertical conductance is 

193 calculated using only the saturated thickness and properties of the 

194 overlying cell if the head in the underlying cell is below its top. 

195 Default is False. 

196 perched: ({True, False}, optional) 

197 keyword to indicate that when a cell is overlying a dewatered 

198 convertible cell, the head difference used in Darcy’s Law is equal to 

199 the head in the overlying cell minus the bottom elevation of the 

200 overlying cell. If not specified, then the default is to use the head 

201 difference between the two cells. 

202 Default is False. 

203 save_specific_discharge: ({True, False}, optional) 

204 keyword to indicate that x, y, and z components of specific discharge 

205 will be calculated at cell centers and written to the cell-by-cell flow 

206 file, which is specified with"budget save file" in Output Control. If 

207 this option is activated, then additional information may be required in 

208 the discretization packages and the GWF Exchange package (if GWF models 

209 are coupled). Specifically, angldegx must be specified in the 

210 connectiondata block of the disu package; angldegx must also be 

211 specified for the GWF Exchange as an auxiliary variable. disu package 

212 has not been implemented yet. 

213 Default is False. 

214 save_saturation: ({True, False}, optional) 

215 keyword to indicate that cell saturation will be written to the budget 

216 file, which is specified with "BUDGET SAVE FILE" in Output Control. 

217 Saturation will be saved to the budget file as an auxiliary variable 

218 saved with the DATA-SAT text label. Saturation is a cell variable that 

219 ranges from zero to one and can be used by post processing programs to 

220 determine how much of a cell volume is saturated. If ICELLTYPE is 0, 

221 then saturation is always one. 

222 xt3d_option: ({True, False}, optional) 

223 If True, the XT3D formulation will be used. By default False. 

224 rhs_option: ({True, False}, optional) 

225 If True, then the XT3D additional terms will be added to the right-hand 

226 side. If False, then the XT3D terms will be put into the coefficient 

227 matrix. By default False. 

228 validate: {True, False} 

229 Flag to indicate whether the package should be validated upon 

230 initialization. This raises a ValidationError if package input is 

231 provided in the wrong manner. Defaults to True. 

232 """ 

233 

234 _pkg_id = "npf" 

235 

236 _init_schemata = { 

237 "icelltype": [ 

238 DTypeSchema(np.integer), 

239 IndexesSchema(), 

240 PKG_DIMS_SCHEMA, 

241 ], 

242 "k": [ 

243 DTypeSchema(np.floating), 

244 IndexesSchema(), 

245 PKG_DIMS_SCHEMA, 

246 ], 

247 "rewet_layer": [ 

248 DTypeSchema(np.floating), 

249 IndexesSchema(), 

250 PKG_DIMS_SCHEMA, 

251 ], 

252 "k22": [ 

253 DTypeSchema(np.floating), 

254 IndexesSchema(), 

255 PKG_DIMS_SCHEMA, 

256 ], 

257 "k33": [ 

258 DTypeSchema(np.floating), 

259 IndexesSchema(), 

260 PKG_DIMS_SCHEMA, 

261 ], 

262 "angle1": [ 

263 DTypeSchema(np.floating), 

264 IndexesSchema(), 

265 PKG_DIMS_SCHEMA, 

266 ], 

267 "angle2": [ 

268 DTypeSchema(np.floating), 

269 IndexesSchema(), 

270 PKG_DIMS_SCHEMA, 

271 ], 

272 "angle3": [ 

273 DTypeSchema(np.floating), 

274 IndexesSchema(), 

275 PKG_DIMS_SCHEMA, 

276 ], 

277 "alternative_cell_averaging": [ 

278 DTypeSchema(str), 

279 DimsSchema(), 

280 CompatibleSettingsSchema("xt3d_option", False), 

281 ], 

282 "rhs_option": [ 

283 DTypeSchema(np.bool_), 

284 DimsSchema(), 

285 CompatibleSettingsSchema("xt3d_option", True), 

286 ], 

287 "save_flows": [DTypeSchema(np.bool_), DimsSchema()], 

288 "starting_head_as_confined_thickness": [DTypeSchema(np.bool_)], 

289 "variable_vertical_conductance": [DTypeSchema(np.bool_)], 

290 "dewatered": [DTypeSchema(np.bool_)], 

291 "perched": [DTypeSchema(np.bool_)], 

292 "save_specific_discharge": [DTypeSchema(np.bool_), DimsSchema()], 

293 } 

294 

295 _write_schemata = { 

296 "k": ( 

297 AllValueSchema(">", 0.0), 

298 IdentityNoDataSchema(other="idomain", is_other_notnull=(">", 0)), 

299 ), 

300 "rewet_layer": ( 

301 IdentityNoDataSchema(other="idomain", is_other_notnull=(">", 0)), 

302 ), 

303 "k22": ( 

304 AllValueSchema(">", 0.0), 

305 IdentityNoDataSchema(other="idomain", is_other_notnull=(">", 0)), 

306 # No need to check coords: dataset ensures they align with idomain. 

307 ), 

308 "k33": ( 

309 AllValueSchema(">", 0.0), 

310 IdentityNoDataSchema(other="idomain", is_other_notnull=(">", 0)), 

311 # No need to check coords: dataset ensures they align with idomain. 

312 ), 

313 "angle1": (IdentityNoDataSchema(other="idomain", is_other_notnull=(">", 0)),), 

314 "angle2": (IdentityNoDataSchema(other="idomain", is_other_notnull=(">", 0)),), 

315 "angle3": (IdentityNoDataSchema(other="idomain", is_other_notnull=(">", 0)),), 

316 } 

317 

318 _grid_data = { 

319 "icelltype": np.int32, 

320 "k": np.float64, 

321 "rewet_layer": np.float64, 

322 "k22": np.float64, 

323 "k33": np.float64, 

324 "angle1": np.float64, 

325 "angle2": np.float64, 

326 "angle3": np.float64, 

327 } 

328 _keyword_map = { 

329 "rewet": "rewet_record", 

330 "rewet_factor": "wetfct", 

331 "rewet_method": "ihdwet", 

332 "rewet_layer": "wetdry", 

333 "variable_vertical_conductance": "variablecv", 

334 "starting_head_as_confined_thickness": "thickstrt", 

335 "rewet_iterations": "iwetit", 

336 "xt3d_option": "xt3doptions", 

337 "rhs_option": "rhs", 

338 } 

339 _template = Package._initialize_template(_pkg_id) 

340 

341 _regrid_method = { 

342 "icelltype": (RegridderType.OVERLAP, "mean"), 

343 "k": (RegridderType.OVERLAP, "geometric_mean"), # horizontal if angle2 = 0 

344 "k22": ( 

345 RegridderType.OVERLAP, 

346 "geometric_mean", 

347 ), # horizontal if angle2 = 0 & angle3 = 0 

348 "k33": ( 

349 RegridderType.OVERLAP, 

350 "harmonic_mean", 

351 ), # vertical if angle2 = 0 & angle3 = 0 

352 "angle1": (RegridderType.OVERLAP, "mean"), 

353 "angle2": (RegridderType.OVERLAP, "mean"), 

354 "angle3": (RegridderType.OVERLAP, "mean"), 

355 "rewet_layer": (RegridderType.OVERLAP, "mean"), 

356 } 

357 

358 @init_log_decorator() 

359 def __init__( 

360 self, 

361 icelltype, 

362 k, 

363 rewet=False, 

364 rewet_layer=None, 

365 rewet_factor=None, 

366 rewet_iterations=None, 

367 rewet_method=None, 

368 k22=None, 

369 k33=None, 

370 angle1=None, 

371 angle2=None, 

372 angle3=None, 

373 cell_averaging=None, 

374 alternative_cell_averaging=None, 

375 save_flows=False, 

376 starting_head_as_confined_thickness=False, 

377 variable_vertical_conductance=False, 

378 dewatered=False, 

379 perched=False, 

380 save_specific_discharge=False, 

381 save_saturation=False, 

382 xt3d_option=False, 

383 rhs_option=False, 

384 validate: bool = True, 

385 ): 

386 # check rewetting 

387 if not rewet and any( 

388 [rewet_layer, rewet_factor, rewet_iterations, rewet_method] 

389 ): 

390 raise ValueError( 

391 "rewet_layer, rewet_factor, rewet_iterations, and rewet_method should" 

392 " all be left at a default value of None if rewet is False." 

393 ) 

394 if cell_averaging is not None: 

395 warnings.warn( 

396 "Use of `cell_averaging` is deprecated, please use `alternative_cell_averaging` instead", 

397 DeprecationWarning, 

398 ) 

399 alternative_cell_averaging = cell_averaging 

400 

401 dict_dataset = { 

402 "icelltype": icelltype, 

403 "k": k, 

404 "rewet": rewet, 

405 "rewet_layer": rewet_layer, 

406 "rewet_factor": rewet_factor, 

407 "rewet_iterations": rewet_iterations, 

408 "rewet_method": rewet_method, 

409 "k22": k22, 

410 "k33": k33, 

411 "angle1": angle1, 

412 "angle2": angle2, 

413 "angle3": angle3, 

414 "alternative_cell_averaging": alternative_cell_averaging, 

415 "save_flows": save_flows, 

416 "starting_head_as_confined_thickness": starting_head_as_confined_thickness, 

417 "variable_vertical_conductance": variable_vertical_conductance, 

418 "dewatered": dewatered, 

419 "perched": perched, 

420 "save_specific_discharge": save_specific_discharge, 

421 "save_saturation": save_saturation, 

422 "xt3d_option": xt3d_option, 

423 "rhs_option": rhs_option, 

424 } 

425 super().__init__(dict_dataset) 

426 self._validate_init_schemata(validate) 

427 

428 def get_xt3d_option(self) -> bool: 

429 """ 

430 Returns the xt3d option value for this object. 

431 """ 

432 return _dataarray_to_bool(self.dataset["xt3d_option"]) 

433 

434 def set_xt3d_option(self, is_xt3d_used: bool, is_rhs: bool) -> None: 

435 """ 

436 Returns the xt3d option value for this object. 

437 """ 

438 self.dataset["rhs_option"] = is_rhs 

439 self.dataset["xt3d_option"] = is_xt3d_used 

440 

441 @property 

442 def is_variable_vertical_conductance(self) -> bool: 

443 """ 

444 Returns the VariableCV option value for this object. 

445 """ 

446 return _dataarray_to_bool(self.dataset["variable_vertical_conductance"]) 

447 

448 @property 

449 def is_dewatered(self) -> bool: 

450 """ 

451 Returns the "dewatered" option value for this object. Used only when variable_vertical_conductance is true 

452 """ 

453 return _dataarray_to_bool(self.dataset["dewatered"]) 

454 

455 def _validate(self, schemata, **kwargs): 

456 # Insert additional kwargs 

457 kwargs["xt3d_option"] = self["xt3d_option"] 

458 errors = super()._validate(schemata, **kwargs) 

459 

460 return errors 

461 

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

463 return self._regrid_method