Coverage for fixtures\mf6_twri_fixture.py: 100%

116 statements  

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

1import numpy as np 

2import pandas as pd 

3import pytest 

4import xarray as xr 

5 

6import imod 

7from imod.typing.grid import zeros_like 

8 

9 

10def make_twri_model(): 

11 nlay = 3 

12 nrow = 15 

13 ncol = 15 

14 shape = (nlay, nrow, ncol) 

15 

16 dx = 5000.0 

17 dy = -5000.0 

18 xmin = 0.0 

19 xmax = dx * ncol 

20 ymin = 0.0 

21 ymax = abs(dy) * nrow 

22 dims = ("layer", "y", "x") 

23 

24 layer = np.array([1, 2, 3]) 

25 y = np.arange(ymax, ymin, dy) + 0.5 * dy 

26 x = np.arange(xmin, xmax, dx) + 0.5 * dx 

27 coords = {"layer": layer, "y": y, "x": x} 

28 

29 # Discretization data 

30 like = xr.DataArray(np.ones(shape), coords=coords, dims=dims) 

31 idomain = like.astype(np.int8) 

32 bottom = xr.DataArray([-200.0, -300.0, -450.0], {"layer": layer}, ("layer",)) 

33 

34 # Constant head 

35 head = xr.full_like(like, np.nan).sel(layer=[1, 2]) 

36 head[..., 0] = 0.0 

37 

38 # Drainage 

39 elevation = xr.full_like(like.sel(layer=1), np.nan) 

40 conductance = xr.full_like(like.sel(layer=1), np.nan) 

41 elevation[7, 1:10] = np.array([0.0, 0.0, 10.0, 20.0, 30.0, 50.0, 70.0, 90.0, 100.0]) 

42 conductance[7, 1:10] = 1.0 

43 

44 # Node properties 

45 icelltype = xr.DataArray([1, 0, 0], {"layer": layer}, ("layer",)) 

46 k = xr.DataArray([1.0e-3, 1.0e-4, 2.0e-4], {"layer": layer}, ("layer",)) 

47 k33 = xr.DataArray([2.0e-8, 2.0e-8, 2.0e-8], {"layer": layer}, ("layer",)) 

48 

49 # Recharge 

50 rch_rate = xr.full_like(like.sel(layer=1), 3.0e-8) 

51 

52 # Well 

53 wells_x = [ 

54 52500.0, 

55 27500.0, 

56 57500.0, 

57 37500.0, 

58 47500.0, 

59 57500.0, 

60 67500.0, 

61 37500.0, 

62 47500.0, 

63 57500.0, 

64 67500.0, 

65 37500.0, 

66 47500.0, 

67 57500.0, 

68 67500.0, 

69 ] 

70 wells_y = [ 

71 52500.0, 

72 57500.0, 

73 47500.0, 

74 32500.0, 

75 32500.0, 

76 32500.0, 

77 32500.0, 

78 22500.0, 

79 22500.0, 

80 22500.0, 

81 22500.0, 

82 12500.0, 

83 12500.0, 

84 12500.0, 

85 12500.0, 

86 ] 

87 screen_top = [ 

88 -300.0, 

89 -200.0, 

90 -200.0, 

91 200.0, 

92 200.0, 

93 200.0, 

94 200.0, 

95 200.0, 

96 200.0, 

97 200.0, 

98 200.0, 

99 200.0, 

100 200.0, 

101 200.0, 

102 200.0, 

103 ] 

104 screen_bottom = [ 

105 -450.0, 

106 -300.0, 

107 -300.0, 

108 -200.0, 

109 -200.0, 

110 -200.0, 

111 -200.0, 

112 -200.0, 

113 -200.0, 

114 -200.0, 

115 -200.0, 

116 -200.0, 

117 -200.0, 

118 -200.0, 

119 -200.0, 

120 ] 

121 rate_wel = [ 

122 -5.0, 

123 -5.0, 

124 -5.0, 

125 -5.0, 

126 -5.0, 

127 -5.0, 

128 -5.0, 

129 -5.0, 

130 -5.0, 

131 -5.0, 

132 -5.0, 

133 -5.0, 

134 -5.0, 

135 -5.0, 

136 -5.0, 

137 ] 

138 

139 # Create and fill the groundwater model. 

140 gwf_model = imod.mf6.GroundwaterFlowModel() 

141 gwf_model["wel"] = imod.mf6.Well( 

142 x=wells_x, 

143 y=wells_y, 

144 screen_top=screen_top, 

145 screen_bottom=screen_bottom, 

146 rate=rate_wel, 

147 minimum_k=1e-19, 

148 save_flows=True, 

149 ) 

150 gwf_model["dis"] = imod.mf6.StructuredDiscretization( 

151 top=200.0, bottom=bottom, idomain=idomain 

152 ) 

153 gwf_model["chd"] = imod.mf6.ConstantHead( 

154 head, print_input=True, print_flows=True, save_flows=True 

155 ) 

156 gwf_model["drn"] = imod.mf6.Drainage( 

157 elevation=elevation, 

158 conductance=conductance, 

159 print_input=True, 

160 print_flows=True, 

161 save_flows=True, 

162 ) 

163 gwf_model["ic"] = imod.mf6.InitialConditions(start=0.0) 

164 gwf_model["npf"] = imod.mf6.NodePropertyFlow( 

165 icelltype=icelltype, 

166 k=k, 

167 k33=k33, 

168 variable_vertical_conductance=True, 

169 dewatered=True, 

170 perched=True, 

171 save_flows=True, 

172 ) 

173 gwf_model["oc"] = imod.mf6.OutputControl(save_head="all", save_budget="all") 

174 gwf_model["rch"] = imod.mf6.Recharge(rch_rate) 

175 

176 gwf_model["sto"] = imod.mf6.SpecificStorage( 

177 specific_storage=1.0e-15, 

178 specific_yield=0.15, 

179 convertible=0, 

180 transient=False, 

181 ) 

182 

183 # The mere presence of the api package should not impact results. It 

184 # should not prevent model execution, splitting or partitioning. 

185 gwf_model["api"] = imod.mf6.ApiPackage(maxbound=20) 

186 

187 # Attach it to a simulation 

188 simulation = imod.mf6.Modflow6Simulation("ex01-twri") 

189 simulation["GWF_1"] = gwf_model 

190 # Define solver settings 

191 simulation["solver"] = imod.mf6.Solution( 

192 modelnames=["GWF_1"], 

193 print_option="summary", 

194 csv_output=False, 

195 no_ptc=True, 

196 outer_dvclose=1.0e-4, 

197 outer_maximum=500, 

198 under_relaxation=None, 

199 inner_dvclose=1.0e-4, 

200 inner_rclose=0.001, 

201 inner_maximum=100, 

202 linear_acceleration="cg", 

203 scaling_method=None, 

204 reordering_method=None, 

205 relaxation_factor=0.97, 

206 ) 

207 # Collect time discretization 

208 simulation.create_time_discretization(additional_times=["2000-01-01", "2000-01-02"]) 

209 return simulation 

210 

211 

212@pytest.fixture(scope="function") 

213def twri_model(): 

214 """Returns steady-state confined model.""" 

215 return make_twri_model() 

216 

217 

218@pytest.fixture(scope="function") 

219def transient_twri_model(): 

220 """Returns transient confined model.""" 

221 simulation = make_twri_model() 

222 gwf_model = simulation["GWF_1"] 

223 like = gwf_model["dis"]["idomain"].astype(float) 

224 gwf_model["sto"] = imod.mf6.SpecificStorage( 

225 specific_storage=xr.full_like(like, 1.0e-15), 

226 specific_yield=xr.full_like(like, 0.15), 

227 convertible=0, 

228 transient=True, 

229 save_flows=True, 

230 ) 

231 simulation.create_time_discretization( 

232 additional_times=pd.date_range("2000-01-01", " 2000-01-31") 

233 ) 

234 return simulation 

235 

236 

237@pytest.fixture(scope="function") 

238def transient_unconfined_twri_model(): 

239 """Returns transient unconfined model, also saves specific discharges.""" 

240 simulation = make_twri_model() 

241 gwf_model = simulation["GWF_1"] 

242 like = gwf_model["dis"]["idomain"].astype(float) 

243 # Force storage to unconfined 

244 gwf_model["sto"] = imod.mf6.SpecificStorage( 

245 specific_storage=xr.full_like(like, 1.0e-15), 

246 specific_yield=xr.full_like(like, 0.15), 

247 convertible=1, 

248 transient=True, 

249 save_flows=True, 

250 ) 

251 # Force npf to unconfined 

252 layer = np.array([1, 2, 3]) 

253 icelltype = xr.DataArray([1, 1, 1], {"layer": layer}, ("layer",)) 

254 gwf_model["npf"]["icelltype"] = icelltype 

255 # Store save cell saturation 

256 gwf_model["npf"]["save_saturation"] = True 

257 # Write specific discharges 

258 gwf_model["npf"]["save_specific_discharge"] = True 

259 simulation.create_time_discretization( 

260 additional_times=pd.date_range("2000-01-01", " 2000-01-31") 

261 ) 

262 return simulation 

263 

264 

265@pytest.mark.usefixtures("twri_model") 

266@pytest.fixture(scope="function") 

267def twri_result(tmpdir_factory): 

268 # Using a tmpdir_factory is the canonical way of sharing a tempory pytest 

269 # directory between different testing modules. 

270 modeldir = tmpdir_factory.mktemp("ex01-twri") 

271 simulation = make_twri_model() 

272 simulation.write(modeldir) 

273 simulation.run() 

274 return modeldir 

275 

276 

277@pytest.mark.usefixtures("transient_twri_model") 

278@pytest.fixture(scope="function") 

279def transient_twri_result(tmpdir_factory, transient_twri_model): 

280 # Using a tmpdir_factory is the canonical way of sharing a tempory pytest 

281 # directory between different testing modules. 

282 modeldir = tmpdir_factory.mktemp("ex01-twri-transient") 

283 simulation = transient_twri_model 

284 simulation.write(modeldir) 

285 simulation.run() 

286 return modeldir 

287 

288 

289@pytest.mark.usefixtures("transient_unconfined_twri_model") 

290@pytest.fixture(scope="function") 

291def transient_unconfined_twri_result(tmpdir_factory, transient_unconfined_twri_model): 

292 # Using a tmpdir_factory is the canonical way of sharing a tempory pytest 

293 # directory between different testing modules. 

294 modeldir = tmpdir_factory.mktemp("ex01-twri-transient-unconfined") 

295 simulation = transient_unconfined_twri_model 

296 simulation.write(modeldir) 

297 simulation.run() 

298 return modeldir 

299 

300 

301@pytest.mark.usefixtures("transient_twri_model") 

302@pytest.fixture(scope="function") 

303def split_transient_twri_model(transient_twri_model): 

304 active = transient_twri_model["GWF_1"].domain.sel(layer=1) 

305 number_partitions = 3 

306 split_location = np.linspace(active.y.min(), active.y.max(), number_partitions + 1) 

307 

308 coords = active.coords 

309 submodel_labels = zeros_like(active) 

310 for id in np.arange(1, number_partitions): 

311 submodel_labels.loc[ 

312 (coords["y"] > split_location[id]) & (coords["y"] <= split_location[id + 1]) 

313 ] = id 

314 

315 split_simulation = transient_twri_model.split(submodel_labels) 

316 

317 return split_simulation