Coverage for C:\src\imod-python\imod\wq\bas.py: 21%

86 statements  

« prev     ^ index     » next       coverage.py v7.4.4, created at 2024-04-08 10:26 +0200

1import jinja2 

2import numpy as np 

3import scipy.ndimage 

4import xarray as xr 

5 

6from imod import util 

7from imod.wq.pkgbase import Package 

8 

9 

10class BasicFlow(Package): 

11 """ 

12 The Basic package is used to specify certain data used in all models. 

13 These include: 

14 

15 1. the locations of acitve, inactive, and specified head in cells, 

16 2. the head stored in inactive cells, 

17 3. the initial head in all cells, and 

18 4. the top and bottom of the aquifer 

19 

20 The number of layers (NLAY) is automatically calculated using the IBOUND. 

21 Thickness is calculated using the specified tops en bottoms. 

22 The Basic package input file is required in all models. 

23 

24 Parameters 

25 ---------- 

26 ibound: xr.DataArray of integers 

27 is the boundary variable. 

28 If IBOUND(J,I,K) < 0, cell J,I,K has a constant head. 

29 If IBOUND(J,I,K) = 0, cell J,I,K is inactive. 

30 If IBOUND(J,I,K) > 0, cell J,I,K is active. 

31 top: float or xr.DataArray of floats 

32 is the top elevation of layer 1. For the common situation in which the 

33 top layer represents a water-table aquifer, it may be reasonable to set 

34 `top` equal to land-surface elevation. 

35 bottom: xr.DataArray of floats 

36 is the bottom elevation of model layers or Quasi-3d confining beds. The 

37 DataArray should at least include the `layer` dimension. 

38 starting_head: float or xr.DataArray of floats 

39 is initial (starting) head—that is, head at the beginning of the 

40 simulation (STRT). starting_head must be specified for all simulations, 

41 including steady-state simulations. One value is read for every model 

42 cell. Usually, these values are read a layer at a time. 

43 inactive_head: float, optional 

44 is the value of head to be assigned to all inactive (no flow) cells 

45 (IBOUND = 0) throughout the simulation (HNOFLO). Because head at 

46 inactive cells is unused in model calculations, this does not affect 

47 model results but serves to identify inactive cells when head is 

48 printed. This value is also used as drawdown at inactive cells if the 

49 drawdown option is used. Even if the user does not anticipate having 

50 inactive cells, a value for inactive_head must be entered. Default 

51 value is 1.0e30. 

52 """ 

53 

54 _pkg_id = "bas6" 

55 _template = jinja2.Template( 

56 "[bas6]\n" 

57 " {%- for layer, value in ibound.items() %}\n" 

58 " ibound_l{{layer}} = {{value}}\n" 

59 " {%- endfor %}\n" 

60 " hnoflo = {{inactive_head}}\n" 

61 " {%- for layer, value in starting_head.items() %}\n" 

62 " strt_l{{layer}} = {{value}}\n" 

63 " {%- endfor -%}" 

64 ) 

65 

66 def __init__(self, ibound, top, bottom, starting_head, inactive_head=1.0e30): 

67 self._check_ibound(ibound) 

68 super().__init__() 

69 self["ibound"] = ibound 

70 self["top"] = top 

71 self["bottom"] = bottom 

72 self["starting_head"] = starting_head 

73 self["inactive_head"] = inactive_head 

74 

75 def _check_ibound(self, ibound): 

76 if not isinstance(ibound, xr.DataArray): 

77 raise TypeError("ibound must be xarray.DataArray") 

78 dims = ibound.dims 

79 if not (dims == ("layer", "y", "x") or dims == ("z", "y", "x")): 

80 raise ValueError( 

81 f'ibound dimensions must be ("layer", "y", "x") or ("z", "y", "x"),' 

82 f" got instead {dims}" 

83 ) 

84 

85 def _render(self, directory, nlayer, *args, **kwargs): 

86 """ 

87 Renders part of runfile that ends up under [bas] section. 

88 """ 

89 d = {} 

90 for varname in ("ibound", "starting_head"): 

91 d[varname] = self._compose_values_layer(varname, directory, nlayer) 

92 d["inactive_head"] = self["inactive_head"].values 

93 

94 return self._template.render(d) 

95 

96 def _compose_top(self, directory): 

97 """ 

98 Composes paths to file, or gets the appropriate scalar value for 

99 a top of model domain. 

100 

101 Parameters 

102 ---------- 

103 directory : str 

104 """ 

105 da = self["top"] 

106 if "x" in da.coords and "y" in da.coords: 

107 if not len(da.shape) == 2: 

108 raise ValueError("Top should either be 2d or a scalar value") 

109 d = {} 

110 d["name"] = "top" 

111 d["directory"] = directory 

112 d["extension"] = ".idf" 

113 value = self._compose_path(d) 

114 else: 

115 if not da.shape == (): 

116 raise ValueError("Top should either be 2d or a scalar value") 

117 value = float(da) 

118 return value 

119 

120 @staticmethod 

121 def _cellsizes(dx): 

122 ncell = dx.size 

123 index_ends = np.argwhere(np.diff(dx) != 0.0) + 1 

124 index_ends = np.append(index_ends, ncell) 

125 index_starts = np.insert(index_ends[:-1], 0, 0) + 1 

126 

127 d = {} 

128 for s, e in zip(index_starts, index_ends): 

129 value = abs(float(dx[s - 1])) 

130 if s == e: 

131 d[f"{s}"] = value 

132 else: 

133 d[f"{s}:{e}"] = value 

134 return d 

135 

136 def _render_dis(self, directory, nlayer): 

137 """ 

138 Renders part of runfile that ends up under [dis] section. 

139 """ 

140 d = {} 

141 d["top"] = self._compose_top(directory) 

142 d["bottom"] = self._compose_values_layer("bottom", directory, nlayer) 

143 d["nlay"], d["nrow"], d["ncol"] = self["ibound"].shape 

144 # TODO: check dx > 0, dy < 0? 

145 if "dx" not in self.dataset or "dy" not in self.dataset: # assume equidistant 

146 dx, _, _ = util.spatial.coord_reference(self.dataset["x"]) 

147 dy, _, _ = util.spatial.coord_reference(self.dataset["y"]) 

148 else: 

149 dx = self.dataset.coords["dx"] 

150 dy = self.dataset.coords["dy"] 

151 

152 if isinstance(dy, (float, int)) or dy.shape in ((), (1,)): 

153 d["dy"] = {"?": abs(float(dy))} 

154 else: 

155 d["dy"] = self._cellsizes(dy) 

156 

157 if isinstance(dx, (float, int)) or dx.shape in ((), (1,)): 

158 d["dx"] = {"?": abs(float(dx))} 

159 else: 

160 d["dx"] = self._cellsizes(dx) 

161 

162 # Non-time dependent part of dis 

163 # Can be inferred from ibound 

164 _dis_template = jinja2.Template( 

165 "[dis]\n" 

166 " nlay = {{nlay}}\n" 

167 " nrow = {{nrow}}\n" 

168 " ncol = {{ncol}}\n" 

169 " {%- for row, value in dy.items() %}\n" 

170 " delc_r{{row}} = {{value}}\n" 

171 " {%- endfor %}\n" 

172 " {%- for col, value in dx.items() %}\n" 

173 " delr_c{{col}} = {{value}}\n" 

174 " {%- endfor %}\n" 

175 " top = {{top}}\n" 

176 " {%- for layer, value in bottom.items() %}\n" 

177 " botm_l{{layer}} = {{value}}\n" 

178 " {%- endfor %}\n" 

179 " laycbd_l? = 0" 

180 ) 

181 

182 return _dis_template.render(d) 

183 

184 def thickness(self): 

185 """ 

186 Computes layer thickness from top and bottom data. 

187 

188 Returns 

189 ------- 

190 thickness : xr.DataArray 

191 """ 

192 th = xr.concat( 

193 [ 

194 self["top"] - self["bottom"].sel(layer=1), 

195 -1.0 * self["bottom"].diff("layer"), 

196 ], 

197 dim="layer", 

198 ) 

199 

200 # Concat results in wrong dimension order (y, x, layer), transpose to 

201 # fix dimension order. 

202 return th.transpose("layer", "y", "x", missing_dims="ignore") 

203 

204 def _pkgcheck(self, ibound=None): 

205 if (self["top"] < self["bottom"]).any(): 

206 raise ValueError(f"top should be larger than bottom in {self}") 

207 

208 active_cells = self["ibound"] != 0 

209 if (active_cells & np.isnan(self["starting_head"])).any(): 

210 raise ValueError( 

211 f"Active cells in ibound may not have a nan value in starting_head in {self}" 

212 ) 

213 

214 _, nlabels = scipy.ndimage.label(active_cells.values) 

215 if nlabels > 1: 

216 raise ValueError( 

217 f"{nlabels} disconnected model domain detected in the ibound in {self}" 

218 )