Coverage for fixtures\mf6_twri_fixture.py: 100%
116 statements
« prev ^ index » next coverage.py v7.5.1, created at 2024-05-13 11:15 +0200
« prev ^ index » next coverage.py v7.5.1, created at 2024-05-13 11:15 +0200
1import numpy as np
2import pandas as pd
3import pytest
4import xarray as xr
6import imod
7from imod.typing.grid import zeros_like
10def make_twri_model():
11 nlay = 3
12 nrow = 15
13 ncol = 15
14 shape = (nlay, nrow, ncol)
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")
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}
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",))
34 # Constant head
35 head = xr.full_like(like, np.nan).sel(layer=[1, 2])
36 head[..., 0] = 0.0
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
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",))
49 # Recharge
50 rch_rate = xr.full_like(like.sel(layer=1), 3.0e-8)
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 ]
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)
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 )
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)
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 outer_dvclose=1.0e-4,
195 outer_maximum=500,
196 under_relaxation=None,
197 inner_dvclose=1.0e-4,
198 inner_rclose=0.001,
199 inner_maximum=100,
200 linear_acceleration="cg",
201 scaling_method=None,
202 reordering_method=None,
203 relaxation_factor=0.97,
204 )
205 # Collect time discretization
206 simulation.create_time_discretization(additional_times=["2000-01-01", "2000-01-02"])
207 return simulation
210@pytest.fixture(scope="function")
211def twri_model():
212 """Returns steady-state confined model."""
213 return make_twri_model()
216@pytest.fixture(scope="function")
217def transient_twri_model():
218 """Returns transient confined model."""
219 simulation = make_twri_model()
220 gwf_model = simulation["GWF_1"]
221 like = gwf_model["dis"]["idomain"].astype(float)
222 gwf_model["sto"] = imod.mf6.SpecificStorage(
223 specific_storage=xr.full_like(like, 1.0e-15),
224 specific_yield=xr.full_like(like, 0.15),
225 convertible=0,
226 transient=True,
227 save_flows=True,
228 )
229 simulation.create_time_discretization(
230 additional_times=pd.date_range("2000-01-01", " 2000-01-31")
231 )
232 return simulation
235@pytest.fixture(scope="function")
236def transient_unconfined_twri_model():
237 """Returns transient unconfined model, also saves specific discharges."""
238 simulation = make_twri_model()
239 gwf_model = simulation["GWF_1"]
240 like = gwf_model["dis"]["idomain"].astype(float)
241 # Force storage to unconfined
242 gwf_model["sto"] = imod.mf6.SpecificStorage(
243 specific_storage=xr.full_like(like, 1.0e-15),
244 specific_yield=xr.full_like(like, 0.15),
245 convertible=1,
246 transient=True,
247 save_flows=True,
248 )
249 # Force npf to unconfined
250 layer = np.array([1, 2, 3])
251 icelltype = xr.DataArray([1, 1, 1], {"layer": layer}, ("layer",))
252 gwf_model["npf"]["icelltype"] = icelltype
253 # Store save cell saturation
254 gwf_model["npf"]["save_saturation"] = True
255 # Write specific discharges
256 gwf_model["npf"]["save_specific_discharge"] = True
257 simulation.create_time_discretization(
258 additional_times=pd.date_range("2000-01-01", " 2000-01-31")
259 )
260 return simulation
263@pytest.mark.usefixtures("twri_model")
264@pytest.fixture(scope="function")
265def twri_result(tmpdir_factory):
266 # Using a tmpdir_factory is the canonical way of sharing a tempory pytest
267 # directory between different testing modules.
268 modeldir = tmpdir_factory.mktemp("ex01-twri")
269 simulation = make_twri_model()
270 simulation.write(modeldir)
271 simulation.run()
272 return modeldir
275@pytest.mark.usefixtures("transient_twri_model")
276@pytest.fixture(scope="function")
277def transient_twri_result(tmpdir_factory, transient_twri_model):
278 # Using a tmpdir_factory is the canonical way of sharing a tempory pytest
279 # directory between different testing modules.
280 modeldir = tmpdir_factory.mktemp("ex01-twri-transient")
281 simulation = transient_twri_model
282 simulation.write(modeldir)
283 simulation.run()
284 return modeldir
287@pytest.mark.usefixtures("transient_unconfined_twri_model")
288@pytest.fixture(scope="function")
289def transient_unconfined_twri_result(tmpdir_factory, transient_unconfined_twri_model):
290 # Using a tmpdir_factory is the canonical way of sharing a tempory pytest
291 # directory between different testing modules.
292 modeldir = tmpdir_factory.mktemp("ex01-twri-transient-unconfined")
293 simulation = transient_unconfined_twri_model
294 simulation.write(modeldir)
295 simulation.run()
296 return modeldir
299@pytest.mark.usefixtures("transient_twri_model")
300@pytest.fixture(scope="function")
301def split_transient_twri_model(transient_twri_model):
302 active = transient_twri_model["GWF_1"].domain.sel(layer=1)
303 number_partitions = 3
304 split_location = np.linspace(active.y.min(), active.y.max(), number_partitions + 1)
306 coords = active.coords
307 submodel_labels = zeros_like(active)
308 for id in np.arange(1, number_partitions):
309 submodel_labels.loc[
310 (coords["y"] > split_location[id]) & (coords["y"] <= split_location[id + 1])
311 ] = id
313 split_simulation = transient_twri_model.split(submodel_labels)
315 return split_simulation