Coverage for fixtures\mf6_circle_fixture.py: 97%
183 statements
« prev ^ index » next coverage.py v7.5.1, created at 2024-05-16 11:41 +0200
« prev ^ index » next coverage.py v7.5.1, created at 2024-05-16 11:41 +0200
1import copy
3import numpy as np
4import pandas as pd
5import pytest
6import xarray as xr
7import xugrid as xu
9import imod
12def make_circle_model():
13 grid = imod.data.circle()
14 nface = grid.n_face
16 nlayer = 2
18 idomain = xu.UgridDataArray(
19 xr.DataArray(
20 np.ones((nlayer, nface), dtype=np.int32),
21 coords={"layer": [1, 2]},
22 dims=["layer", grid.face_dimension],
23 ),
24 grid=grid,
25 )
26 icelltype = xu.full_like(idomain, 0)
27 k = xu.full_like(idomain, 1.0, dtype=np.float64)
28 k33 = k.copy()
29 rch_rate = xu.full_like(k.sel(layer=1), 0.001, dtype=float)
31 bottom = k * xr.DataArray([5.0, 0.0], dims=["layer"])
32 chd_location = xu.zeros_like(k.sel(layer=2), dtype=bool).ugrid.binary_dilation(
33 border_value=True
34 )
35 constant_head = xu.full_like(k.sel(layer=2), 1.0).where(chd_location)
37 gwf_model = imod.mf6.GroundwaterFlowModel()
38 gwf_model["disv"] = imod.mf6.VerticesDiscretization(
39 top=10.0, bottom=bottom, idomain=idomain
40 )
41 gwf_model["chd"] = imod.mf6.ConstantHead(
42 constant_head, print_input=True, print_flows=True, save_flows=True
43 )
44 gwf_model["ic"] = imod.mf6.InitialConditions(start=0.0)
45 gwf_model["npf"] = imod.mf6.NodePropertyFlow(
46 icelltype=icelltype,
47 k=k,
48 k33=k33,
49 save_flows=True,
50 )
51 gwf_model["sto"] = imod.mf6.SpecificStorage(
52 specific_storage=1.0e-5,
53 specific_yield=0.15,
54 transient=False,
55 convertible=0,
56 save_flows=False,
57 )
58 gwf_model["oc"] = imod.mf6.OutputControl(save_head="last", save_budget="last")
59 gwf_model["rch"] = imod.mf6.Recharge(rch_rate)
61 simulation = imod.mf6.Modflow6Simulation("circle")
62 simulation["GWF_1"] = gwf_model
63 simulation["solver"] = imod.mf6.Solution(
64 modelnames=["GWF_1"],
65 print_option="summary",
66 outer_dvclose=1.0e-4,
67 outer_maximum=500,
68 under_relaxation=None,
69 inner_dvclose=1.0e-4,
70 inner_rclose=0.001,
71 inner_maximum=100,
72 linear_acceleration="cg",
73 scaling_method=None,
74 reordering_method=None,
75 relaxation_factor=0.97,
76 )
77 simtimes = pd.date_range(start="2000-01-01", end="2000-01-03")
78 simulation.create_time_discretization(additional_times=simtimes)
79 return simulation
82def make_circle_model_flow_with_transport_data(species: list[str]):
83 grid = imod.data.circle()
84 max_concentration = 35.0
85 min_concentration = 0.0
86 nface = grid.n_face
88 nlayer = 2
90 idomain = xu.UgridDataArray(
91 xr.DataArray(
92 np.ones((nlayer, nface), dtype=np.int32),
93 coords={"layer": [1, 2]},
94 dims=["layer", grid.face_dimension],
95 ),
96 grid=grid,
97 )
98 icelltype = xu.full_like(idomain, 0)
99 k = xu.full_like(idomain, 1.0, dtype=np.float64)
100 k33 = k.copy()
101 rch_rate = xu.full_like(k.sel(layer=1), 0.001, dtype=float)
102 rch_concentration = xu.full_like(rch_rate, min_concentration)
103 rch_concentration = rch_concentration.expand_dims(species=species)
104 bottom = k * xr.DataArray([5.0, 0.0], dims=["layer"])
105 chd_location = xu.zeros_like(k.sel(layer=2), dtype=bool).ugrid.binary_dilation(
106 border_value=True
107 )
108 constant_head = xu.full_like(k.sel(layer=2), 1.0).where(chd_location)
109 constant_concentration = xu.full_like(constant_head, max_concentration).where(
110 chd_location
111 )
112 constant_concentration = constant_concentration.expand_dims(species=species)
114 gwf_model = imod.mf6.GroundwaterFlowModel(save_flows=True)
115 gwf_model["disv"] = imod.mf6.VerticesDiscretization(
116 top=10.0, bottom=bottom, idomain=idomain
117 )
118 gwf_model["chd"] = imod.mf6.ConstantHead(
119 constant_head,
120 concentration=constant_concentration,
121 print_input=True,
122 print_flows=True,
123 save_flows=True,
124 )
125 gwf_model["ic"] = imod.mf6.InitialConditions(start=0.0)
126 gwf_model["npf"] = imod.mf6.NodePropertyFlow(
127 icelltype=icelltype,
128 k=k,
129 k33=k33,
130 save_flows=True,
131 )
132 gwf_model["sto"] = imod.mf6.SpecificStorage(
133 specific_storage=1.0e-5,
134 specific_yield=0.15,
135 transient=False,
136 convertible=0,
137 save_flows=False,
138 )
139 gwf_model["oc"] = imod.mf6.OutputControl(save_head="last", save_budget="last")
140 gwf_model["rch"] = imod.mf6.Recharge(
141 rch_rate, save_flows=True, concentration=rch_concentration
142 )
144 simulation = imod.mf6.Modflow6Simulation("circle")
145 simulation["GWF_1"] = gwf_model
146 simulation["solver"] = imod.mf6.Solution(
147 modelnames=["GWF_1"],
148 print_option="summary",
149 outer_dvclose=1.0e-4,
150 outer_maximum=500,
151 under_relaxation=None,
152 inner_dvclose=1.0e-4,
153 inner_rclose=0.001,
154 inner_maximum=100,
155 linear_acceleration="cg",
156 scaling_method=None,
157 reordering_method=None,
158 relaxation_factor=0.97,
159 )
160 simtimes = pd.date_range(start="2000-01-01", end="2000-01-03")
161 simulation.create_time_discretization(simtimes)
162 return simulation
165@pytest.fixture(scope="function")
166def circle_model():
167 return make_circle_model()
170@pytest.mark.usefixtures("circle_model")
171@pytest.fixture(scope="session")
172def circle_result(tmpdir_factory):
173 # Using a tmpdir_factory is the canonical way of sharing a tempory pytest
174 # directory between different testing modules.
175 modeldir = tmpdir_factory.mktemp("circle")
176 simulation = make_circle_model()
177 simulation.write(modeldir)
178 simulation.run()
179 return modeldir
182def make_circle_model_evt():
183 simulation = make_circle_model()
184 gwf_model = simulation["GWF_1"]
186 idomain = gwf_model["disv"].dataset["idomain"]
187 like = idomain.sel(layer=1).astype(np.float64)
188 face_dim = idomain.ugrid.grid.face_dimension
190 rate = xu.full_like(like, 0.001)
191 # Lay surface on chd level
192 surface = xu.full_like(like, 1.0)
193 depth = xu.full_like(like, 2.0)
195 segments = xr.DataArray(
196 data=[1, 2, 3], coords={"segment": [1, 2, 3]}, dims=("segment",)
197 )
198 segments_reversed = segments.copy()
199 segments_reversed.values = [3, 2, 1]
201 proportion_depth = xu.full_like(like, 0.3) * segments
202 proportion_rate = xu.full_like(like, 0.3) * segments_reversed
204 proportion_depth = proportion_depth.transpose("segment", face_dim)
205 proportion_rate = proportion_rate.transpose("segment", face_dim)
207 gwf_model["evt"] = imod.mf6.Evapotranspiration(
208 surface, rate, depth, proportion_rate, proportion_depth
209 )
211 simulation["GWF_1"] = gwf_model
213 return simulation
216@pytest.fixture(scope="session")
217def circle_model_evt():
218 return make_circle_model_evt()
221@pytest.mark.usefixtures("circle_model_evt")
222@pytest.fixture(scope="session")
223def circle_result_evt(tmpdir_factory):
224 # Using a tmpdir_factory is the canonical way of sharing a tempory pytest
225 # directory between different testing modules.
226 modeldir = tmpdir_factory.mktemp("circle_evt")
227 simulation = make_circle_model_evt()
228 simulation.write(modeldir)
229 simulation.run()
230 return modeldir
233def make_circle_model_save_sto():
234 simulation = make_circle_model()
235 gwf_model = simulation["GWF_1"]
237 gwf_model["sto"].dataset["save_flows"] = True
238 return simulation
241@pytest.fixture(scope="session")
242def circle_result_sto(tmpdir_factory):
243 """
244 Circle result with storage fluxes, which are saved as METH1 instead of METH6
245 """
246 # Using a tmpdir_factory is the canonical way of sharing a tempory pytest
247 # directory between different testing modules.
248 modeldir = tmpdir_factory.mktemp("circle_sto")
249 simulation = make_circle_model_save_sto()
250 simulation.write(modeldir)
251 simulation.run()
252 return modeldir
255@pytest.mark.usefixtures("circle_model_evt")
256@pytest.fixture(scope="function")
257def circle_partitioned():
258 simulation = make_circle_model_evt()
260 idomain = simulation["GWF_1"]["disv"].dataset["idomain"]
261 submodel_labels = copy.deepcopy(idomain.sel({"layer": 1}))
263 submodel_labels.values[:67] = 0
264 submodel_labels.values[67:118] = 1
265 submodel_labels.values[118:] = 2
267 return simulation.split(submodel_labels)
270@pytest.fixture(scope="function")
271def circle_model_transport():
272 al = 0.001
273 porosity = 0.3
274 max_concentration = 35.0
275 min_concentration = 0.0
276 max_density = 1025.0
277 min_density = 1000.0
279 simulation = make_circle_model_flow_with_transport_data(["salinity"])
280 gwf_model = simulation["GWF_1"]
282 slope = (max_density - min_density) / (max_concentration - min_concentration)
283 gwf_model["buoyancy"] = imod.mf6.Buoyancy(
284 reference_density=min_density,
285 modelname=["transport"],
286 reference_concentration=[min_concentration],
287 density_concentration_slope=[slope],
288 species=["salinity"],
289 )
290 transport_model = imod.mf6.GroundwaterTransportModel(save_flows=True)
291 transport_model["ssm"] = imod.mf6.SourceSinkMixing.from_flow_model(
292 gwf_model, "salinity", save_flows=True
293 )
294 transport_model["disv"] = gwf_model["disv"]
296 # %%
297 # Now we define some transport packages for simulating the physical processes
298 # of advection, mechanical dispersion, and molecular diffusion dispersion. This
299 # example is transient, and the volume available for storage is the porosity,
300 # in this case 0.10.
302 transport_model["dsp"] = imod.mf6.Dispersion(
303 diffusion_coefficient=1e-4,
304 longitudinal_horizontal=al,
305 transversal_horizontal1=al * 0.1,
306 transversal_vertical=al * 0.01,
307 xt3d_off=False,
308 xt3d_rhs=False,
309 )
310 transport_model["adv"] = imod.mf6.AdvectionUpstream()
311 transport_model["mst"] = imod.mf6.MobileStorageTransfer(porosity, save_flows=True)
313 # %%
314 # Define the maximum concentration as the initial conditions, also output
315 # options for the transport model, and assign the transport model to the
316 # simulation as well.
317 max_concentration = 35.0
318 min_concentration = 0.0
319 transport_model["ic"] = imod.mf6.InitialConditions(start=max_concentration)
320 transport_model["oc"] = imod.mf6.OutputControl(
321 save_concentration="last", save_budget="last"
322 )
324 simulation["transport"] = transport_model
325 simulation["transport_solver"] = imod.mf6.Solution(
326 modelnames=["transport"],
327 print_option="summary",
328 outer_dvclose=1.0e-4,
329 outer_maximum=500,
330 under_relaxation=None,
331 inner_dvclose=1.0e-4,
332 inner_rclose=0.001,
333 inner_maximum=100,
334 linear_acceleration="bicgstab",
335 scaling_method=None,
336 reordering_method=None,
337 relaxation_factor=0.97,
338 )
339 simtimes = pd.date_range(start="2000-01-01", end="2000-01-03")
340 simulation.create_time_discretization(additional_times=simtimes)
341 return simulation
344@pytest.fixture(scope="function")
345def circle_model_transport_multispecies_variable_density():
346 al = 0.001
347 porosity = 0.3
348 max_concentration = 35.0
349 min_concentration = 0.0
350 max_density = 1025.0
351 min_density = 1000.0
352 species = ["salt", "temp"]
354 simulation = make_circle_model_flow_with_transport_data(species)
355 gwf_model = simulation["GWF_1"]
357 for specie in species:
358 transport_model = imod.mf6.GroundwaterTransportModel(save_flows=True)
359 transport_model["ssm"] = imod.mf6.SourceSinkMixing.from_flow_model(
360 gwf_model, specie, save_flows=True
361 )
362 transport_model["disv"] = gwf_model["disv"]
364 # %%
365 # Now we define some transport packages for simulating the physical processes
366 # of advection, mechanical dispersion, and molecular diffusion dispersion. This
367 # example is transient, and the volume available for storage is the porosity,
368 # in this case 0.10.
370 transport_model["dsp"] = imod.mf6.Dispersion(
371 diffusion_coefficient=1e-4,
372 longitudinal_horizontal=al,
373 transversal_horizontal1=al * 0.1,
374 transversal_vertical=al * 0.01,
375 xt3d_off=False,
376 xt3d_rhs=False,
377 )
378 transport_model["adv"] = imod.mf6.AdvectionUpstream()
379 transport_model["mst"] = imod.mf6.MobileStorageTransfer(
380 porosity, save_flows=True
381 )
383 # %% Define the maximum concentration as the initial conditions, also
384 # output options for the transport model, and assign the transport model
385 # to the simulation as well.
386 transport_model["ic"] = imod.mf6.InitialConditions(start=max_concentration)
387 transport_model["oc"] = imod.mf6.OutputControl(
388 save_concentration="all", save_budget="all"
389 )
391 simulation[f"tpt_{specie}"] = transport_model
392 slope = (max_density - min_density) / (max_concentration - min_concentration)
393 modelnames = [f"tpt_{specie}" for specie in species]
394 gwf_model["buoyancy"] = imod.mf6.Buoyancy(
395 reference_density=min_density,
396 modelname=modelnames,
397 reference_concentration=[min_concentration, min_concentration],
398 density_concentration_slope=[slope, slope],
399 species=species,
400 )
402 simulation["transport_solver"] = imod.mf6.Solution(
403 modelnames=modelnames,
404 print_option="summary",
405 outer_dvclose=1.0e-4,
406 outer_maximum=500,
407 under_relaxation=None,
408 inner_dvclose=1.0e-4,
409 inner_rclose=0.001,
410 inner_maximum=100,
411 linear_acceleration="bicgstab",
412 scaling_method=None,
413 reordering_method=None,
414 relaxation_factor=0.97,
415 )
416 simtimes = pd.date_range(start="2000-01-01", end="2000-01-03")
417 simulation.create_time_discretization(additional_times=simtimes)
418 return simulation