Coverage for C:\src\imod-python\imod\mf6\sto.py: 100%

79 statements  

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

1import abc 

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 DimsSchema, 

14 DTypeSchema, 

15 IdentityNoDataSchema, 

16 IndexesSchema, 

17) 

18 

19 

20class Storage(Package): 

21 _pkg_id = "sto_deprecated" 

22 

23 def __init__(*args, **kwargs): 

24 raise NotImplementedError( 

25 r"Storage package has been removed. Use SpecificStorage or StorageCoefficient instead." 

26 ) 

27 

28 

29class StorageBase(Package, IRegridPackage, abc.ABC): 

30 def get_options(self, d): 

31 # Skip both variables in grid_data and "transient". 

32 not_options = list(self._grid_data.keys()) 

33 not_options += "transient" 

34 

35 for varname in self.dataset.data_vars.keys(): # pylint:disable=no-member 

36 if varname in not_options: 

37 continue 

38 v = self.dataset[varname].values[()] 

39 if self._valid(v): # skip None and False 

40 d[varname] = v 

41 return d 

42 

43 def _render_dict(self, directory, pkgname, globaltimes, binary): 

44 d = {} 

45 stodirectory = directory / pkgname 

46 for varname in self._grid_data: 

47 key = self._keyword_map.get(varname, varname) 

48 layered, value = self._compose_values( 

49 self[varname], stodirectory, key, binary=binary 

50 ) 

51 if self._valid(value): # skip False or None 

52 d[f"{key}_layered"], d[key] = layered, value 

53 

54 periods = {} 

55 if "time" in self.dataset["transient"].coords: 

56 package_times = self.dataset["transient"].coords["time"].values 

57 starts = np.searchsorted(globaltimes, package_times) + 1 

58 for i, s in enumerate(starts): 

59 periods[s] = self.dataset["transient"].isel(time=i).values[()] 

60 else: 

61 periods[1] = self.dataset["transient"].values[()] 

62 

63 d["periods"] = periods 

64 

65 d = self.get_options(d) 

66 

67 return d 

68 

69 

70class SpecificStorage(StorageBase): 

71 """ 

72 Storage Package with specific storage. 

73 

74 From wikipedia (https://en.wikipedia.org/wiki/Specific_storage): 

75 

76 "The specific storage is the amount of water that a portion of an aquifer 

77 releases from storage, per unit mass or volume of aquifer, per unit change 

78 in hydraulic head, while remaining fully saturated." 

79 

80 If the STO Package is not included for a model, then storage changes will 

81 not be calculated, and thus, the model will be steady state. Only one STO 

82 Package can be specified for a GWF model. 

83 

84 Parameters 

85 ---------- 

86 specific_storage: array of floats (xr.DataArray) 

87 Is specific storage. Specific storage values must be greater than 

88 or equal to 0. (ss) 

89 specific_yield: array of floats (xr.DataArray) 

90 Is specific yield. Specific yield values must be greater than or 

91 equal to 0. Specific yield does not have to be specified if there are no 

92 convertible cells (convertible=0 in every cell). (sy) 

93 transient: ({True, False}), or a DataArray with a time coordinate and dtype Bool 

94 Boolean to indicate if the model is transient or steady-state. 

95 convertible: array of int (xr.DataArray) 

96 Is a flag for each cell that specifies whether or not a cell is 

97 convertible for the storage calculation. 0 indicates confined storage is 

98 used. >0 indicates confined storage is used when head is above cell top 

99 and a mixed formulation of unconfined and confined storage is used when 

100 head is below cell top. (iconvert) 

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

102 Indicates that storage flow terms will be written to the file specified 

103 with "BUDGET FILEOUT" in Output Control. Default is False. 

104 validate: {True, False} 

105 Flag to indicate whether the package should be validated upon 

106 initialization. This raises a ValidationError if package input is 

107 provided in the wrong manner. Defaults to True. 

108 """ 

109 

110 _pkg_id = "sto" 

111 _grid_data = { 

112 "convertible": np.int32, 

113 "specific_storage": np.float64, 

114 "specific_yield": np.float64, 

115 } 

116 _keyword_map = { 

117 "specific_storage": "ss", 

118 "specific_yield": "sy", 

119 "convertible": "iconvert", 

120 } 

121 

122 _init_schemata = { 

123 "specific_storage": ( 

124 DTypeSchema(np.floating), 

125 IndexesSchema(), 

126 PKG_DIMS_SCHEMA, 

127 ), 

128 "specific_yield": ( 

129 DTypeSchema(np.floating), 

130 IndexesSchema(), 

131 PKG_DIMS_SCHEMA, 

132 ), 

133 "transient": ( 

134 DTypeSchema(np.bool_), 

135 IndexesSchema(), 

136 DimsSchema("time") | DimsSchema(), 

137 ), 

138 "convertible": ( 

139 DTypeSchema(np.integer), 

140 IndexesSchema(), 

141 PKG_DIMS_SCHEMA, 

142 ), 

143 "save_flows": (DTypeSchema(np.bool_), DimsSchema()), 

144 } 

145 

146 _write_schemata = { 

147 "specific_storage": ( 

148 AllValueSchema(">=", 0.0), 

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

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

151 ), 

152 "specific_yield": ( 

153 AllValueSchema(">=", 0.0), 

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

155 ), 

156 "convertible": ( 

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

158 ), 

159 } 

160 

161 _regrid_method = { 

162 "convertible": (RegridderType.OVERLAP, "mode"), 

163 "specific_storage": (RegridderType.OVERLAP, "mean"), 

164 "specific_yield": (RegridderType.OVERLAP, "mean"), 

165 } 

166 

167 _template = Package._initialize_template(_pkg_id) 

168 

169 @init_log_decorator() 

170 def __init__( 

171 self, 

172 specific_storage, 

173 specific_yield, 

174 transient, 

175 convertible, 

176 save_flows: bool = False, 

177 validate: bool = True, 

178 ): 

179 dict_dataset = { 

180 "specific_storage": specific_storage, 

181 "specific_yield": specific_yield, 

182 "convertible": convertible, 

183 "transient": transient, 

184 "save_flows": save_flows, 

185 } 

186 super().__init__(dict_dataset) 

187 self._validate_init_schemata(validate) 

188 

189 def render(self, directory, pkgname, globaltimes, binary): 

190 d = self._render_dict(directory, pkgname, globaltimes, binary) 

191 return self._template.render(d) 

192 

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

194 return self._regrid_method 

195 

196 

197class StorageCoefficient(StorageBase): 

198 """ 

199 Storage Package with a storage coefficient. Be careful, 

200 this is not the same as the specific storage. 

201 

202 From wikipedia (https://en.wikipedia.org/wiki/Specific_storage): 

203 

204 "Storativity or the storage coefficient is the volume of water released 

205 from storage per unit decline in hydraulic head in the aquifer, per 

206 unit area of the aquifer. Storativity is a dimensionless quantity, and 

207 is always greater than 0. 

208 

209 Under confined conditions: 

210 

211 S = Ss * b, where S is the storage coefficient, 

212 Ss the specific storage, and b the aquifer thickness. 

213 

214 Under unconfined conditions: 

215 

216 S = Sy, where Sy is the specific yield" 

217 

218 If the STO Package is not included for a model, then storage changes will 

219 not be calculated, and thus, the model will be steady state. Only one STO 

220 Package can be specified for a GWF model. 

221 

222 Parameters 

223 ---------- 

224 storage_coefficient: array of floats (xr.DataArray) 

225 Is storage coefficient. Storage coefficient values must be greater than 

226 or equal to 0. (ss) 

227 specific_yield: array of floats (xr.DataArray) 

228 Is specific yield. Specific yield values must be greater than or 

229 equal to 0. Specific yield does not have to be specified if there are no 

230 convertible cells (convertible=0 in every cell). (sy) 

231 transient: ({True, False}) 

232 Boolean to indicate if the model is transient or steady-state. 

233 convertible: array of int (xr.DataArray) 

234 Is a flag for each cell that specifies whether or not a cell is 

235 convertible for the storage calculation. 0 indicates confined storage is 

236 used. >0 indicates confined storage is used when head is above cell top 

237 and a mixed formulation of unconfined and confined storage is used when 

238 head is below cell top. (iconvert) 

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

240 Indicates that storage flow terms will be written to the file specified 

241 with "BUDGET FILEOUT" in Output Control. Default is False. 

242 validate: {True, False} 

243 Flag to indicate whether the package should be validated upon 

244 initialization. This raises a ValidationError if package input is 

245 provided in the wrong manner. Defaults to True. 

246 """ 

247 

248 _pkg_id = "sto" 

249 _grid_data = { 

250 "convertible": np.int32, 

251 "storage_coefficient": np.float64, 

252 "specific_yield": np.float64, 

253 } 

254 _keyword_map = { 

255 "storage_coefficient": "ss", 

256 "specific_yield": "sy", 

257 "convertible": "iconvert", 

258 } 

259 

260 _init_schemata = { 

261 "storage_coefficient": ( 

262 DTypeSchema(np.floating), 

263 IndexesSchema(), 

264 PKG_DIMS_SCHEMA, 

265 ), 

266 "specific_yield": ( 

267 DTypeSchema(np.floating), 

268 IndexesSchema(), 

269 PKG_DIMS_SCHEMA, 

270 ), 

271 "transient": ( 

272 DTypeSchema(np.bool_), 

273 IndexesSchema(), 

274 DimsSchema("time") | DimsSchema(), 

275 ), 

276 "convertible": ( 

277 DTypeSchema(np.integer), 

278 IndexesSchema(), 

279 PKG_DIMS_SCHEMA, 

280 ), 

281 "save_flows": (DTypeSchema(np.bool_), DimsSchema()), 

282 } 

283 

284 _write_schemata = { 

285 "storage_coefficient": ( 

286 AllValueSchema(">=", 0.0), 

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

288 ), 

289 "specific_yield": ( 

290 AllValueSchema(">=", 0.0), 

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

292 ), 

293 "convertible": ( 

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

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

296 ), 

297 } 

298 

299 _regrid_method = { 

300 "convertible": (RegridderType.OVERLAP, "mode"), 

301 "storage_coefficient": (RegridderType.OVERLAP, "mean"), 

302 "specific_yield": (RegridderType.OVERLAP, "mean"), 

303 } 

304 

305 _template = Package._initialize_template(_pkg_id) 

306 

307 @init_log_decorator() 

308 def __init__( 

309 self, 

310 storage_coefficient, 

311 specific_yield, 

312 transient, 

313 convertible, 

314 save_flows: bool = False, 

315 validate: bool = True, 

316 ): 

317 dict_dataset = { 

318 "storage_coefficient": storage_coefficient, 

319 "specific_yield": specific_yield, 

320 "convertible": convertible, 

321 "transient": transient, 

322 "save_flows": save_flows, 

323 } 

324 super().__init__(dict_dataset) 

325 self._validate_init_schemata(validate) 

326 

327 def render(self, directory, pkgname, globaltimes, binary): 

328 d = self._render_dict(directory, pkgname, globaltimes, binary) 

329 d["storagecoefficient"] = True 

330 return self._template.render(d) 

331 

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

333 return self._regrid_method