Coverage for fixtures\package_instance_creation.py: 100%
70 statements
« prev ^ index » next coverage.py v7.4.4, created at 2024-04-08 10:26 +0200
« prev ^ index » next coverage.py v7.4.4, created at 2024-04-08 10:26 +0200
1import geopandas as gpd
2import numpy as np
3import pandas as pd
4import shapely
5import xarray as xr
6import xugrid as xu
8import imod
9import imod.tests.fixtures.mf6_lake_package_fixture as mf_lake
10from imod.mf6 import GWFGWF
12"""
13This file is used to create instances of imod packages for testing purposes.
14The main usage is importing ALL_PACKAGE_INSTANCES into a test- this list contains an instance of
15each packages and boundary condition in mf6.
16"""
19def get_structured_grid_da(dtype, value=1):
20 """
21 This function creates a dataarray with scalar values for a grid of 3 layers and 9 rows and columns.
22 """
23 shape = nlay, nrow, ncol = 3, 9, 9
24 dims = ("layer", "y", "x")
26 dx = 10.0
27 dy = -10.0
28 xmin = 0.0
29 xmax = dx * ncol
30 ymin = 0.0
31 ymax = abs(dy) * nrow
33 layer = np.arange(1, nlay + 1)
34 y = np.arange(ymax, ymin, dy) + 0.5 * dy
35 x = np.arange(xmin, xmax, dx) + 0.5 * dx
36 coords = {"layer": layer, "y": y, "x": x}
38 da = xr.DataArray(np.ones(shape, dtype=dtype) * value, coords=coords, dims=dims)
39 return da
42def get_unstructured_grid_da(dtype, value=1):
43 """
44 This function creates an xugrid dataarray with scalar values for an unstructured grid
45 """
46 grid = imod.data.circle()
47 nface = grid.n_face
48 nlayer = 2
50 dims = ("layer", grid.face_dimension)
51 shape = (nlayer, nface)
53 uda = xu.UgridDataArray(
54 xr.DataArray(
55 np.ones(shape, dtype=dtype) * value,
56 coords={"layer": [1, 2]},
57 dims=dims,
58 ),
59 grid=grid,
60 )
61 return uda
64def get_grid_da(is_unstructured, dtype, value=1):
65 """
66 helper function for creating an xarray dataset of a given type
67 Depending on the is_unstructured input parameter, it will create an array for a
68 structured grid or for an unstructured grid.
69 """
71 if is_unstructured:
72 return get_unstructured_grid_da(dtype, value)
73 else:
74 return get_structured_grid_da(dtype, value)
77def create_lake_package(is_unstructured):
78 is_lake1 = get_grid_da(is_unstructured, bool, False)
79 times_of_numeric_timeseries = [
80 np.datetime64("2000-01-01"),
81 np.datetime64("2000-03-01"),
82 np.datetime64("2000-05-01"),
83 ]
84 numeric = xr.DataArray(
85 np.full((len(times_of_numeric_timeseries)), 5.0),
86 coords={"time": times_of_numeric_timeseries},
87 dims=["time"],
88 )
89 if is_unstructured:
90 is_lake1[1, 0] = True
91 is_lake1[1, 1] = True
92 lake1 = mf_lake.create_lake_data_unstructured(
93 is_lake1, 11.0, "Naardermeer", stage=numeric
94 )
95 else:
96 is_lake1[1, 1, 0] = True
97 is_lake1[1, 1, 1] = True
98 lake1 = mf_lake.create_lake_data_structured(
99 is_lake1, 11.0, "Naardermeer", stage=numeric
100 )
101 outlet1 = imod.mf6.OutletManning("Naardermeer", "", 3.0, 2.0, 3.0, 4.0)
102 return imod.mf6.Lake.from_lakes_and_outlets([lake1], [outlet1])
105def create_vertices_discretization():
106 """
107 return imod.mf6.VerticesDiscretization object
108 """
109 idomain = get_grid_da(True, int, value=1)
110 bottom = idomain * xr.DataArray([5.0, 0.0], dims=["layer"])
111 return imod.mf6.VerticesDiscretization(top=10.0, bottom=bottom, idomain=idomain)
114def create_instance_packages(is_unstructured):
115 """
116 creates instances of those modflow packages that are not boundary conditions.
117 """
118 return [
119 imod.mf6.Dispersion(
120 diffusion_coefficient=get_grid_da(is_unstructured, np.float32, 1e-4),
121 longitudinal_horizontal=get_grid_da(is_unstructured, np.float32, 10),
122 transversal_horizontal1=get_grid_da(is_unstructured, np.float32, 10),
123 longitudinal_vertical=get_grid_da(is_unstructured, np.float32, 5),
124 transversal_horizontal2=get_grid_da(is_unstructured, np.float32, 2),
125 transversal_vertical=get_grid_da(is_unstructured, np.float32, 4),
126 ),
127 imod.mf6.InitialConditions(start=get_grid_da(is_unstructured, np.float32)),
128 imod.mf6.MobileStorageTransfer(
129 porosity=get_grid_da(is_unstructured, np.float32, 0.35),
130 decay=get_grid_da(is_unstructured, np.float32, 0.01),
131 decay_sorbed=get_grid_da(is_unstructured, np.float32, 0.02),
132 bulk_density=get_grid_da(is_unstructured, np.float32, 1300),
133 distcoef=get_grid_da(is_unstructured, np.float32, 0.1),
134 ),
135 imod.mf6.NodePropertyFlow(
136 get_grid_da(is_unstructured, np.int32), 3.0, True, 32.0, 34.0, 7
137 ),
138 imod.mf6.SpecificStorage(
139 0.001, 0.1, True, get_grid_da(is_unstructured, np.int32)
140 ),
141 imod.mf6.StorageCoefficient(
142 0.001, 0.1, True, get_grid_da(is_unstructured, np.int32)
143 ),
144 imod.mf6.ApiPackage(
145 maxbound=33, print_input=True, print_flows=True, save_flows=True
146 ),
147 ]
150def create_instance_boundary_condition_packages(is_unstructured):
151 """
152 creates instances of those modflow packages that are boundary conditions.
153 """
154 return [
155 imod.mf6.ConstantConcentration(
156 get_grid_da(is_unstructured, np.float32, 2),
157 print_input=True,
158 print_flows=True,
159 save_flows=True,
160 ),
161 imod.mf6.ConstantHead(
162 get_grid_da(is_unstructured, np.float32, 2),
163 print_input=True,
164 print_flows=True,
165 save_flows=True,
166 ),
167 imod.mf6.Drainage(
168 elevation=get_grid_da(is_unstructured, np.float64, 4),
169 conductance=get_grid_da(is_unstructured, np.float64, 1e-3),
170 ),
171 imod.mf6.Evapotranspiration(
172 surface=get_grid_da(is_unstructured, np.float64, 3),
173 rate=get_grid_da(is_unstructured, np.float64, 2),
174 depth=get_grid_da(is_unstructured, np.float64, 1),
175 proportion_rate=get_grid_da(is_unstructured, np.float64, 0.2),
176 proportion_depth=get_grid_da(is_unstructured, np.float64, 0.2),
177 fixed_cell=True,
178 ),
179 imod.mf6.GeneralHeadBoundary(
180 head=get_grid_da(is_unstructured, np.float64, 3),
181 conductance=get_grid_da(is_unstructured, np.float64, 0.33),
182 ),
183 imod.mf6.Recharge(
184 rate=get_grid_da(is_unstructured, np.float64, 0.33),
185 ),
186 imod.mf6.River(
187 stage=get_grid_da(is_unstructured, np.float64, 0.33),
188 conductance=get_grid_da(is_unstructured, np.float64, 0.33),
189 bottom_elevation=get_grid_da(is_unstructured, np.float64, 0.33),
190 ),
191 imod.mf6.SourceSinkMixing(
192 package_names=["a", "b"],
193 concentration_boundary_type=["a", "b"],
194 auxiliary_variable_name=["a", "b"],
195 ),
196 create_lake_package(is_unstructured),
197 ]
200STRUCTURED_GRID_PACKAGES = [
201 imod.mf6.StructuredDiscretization(
202 2.0, get_grid_da(False, np.float32), get_grid_da(False, np.int32)
203 ),
204 imod.mf6.WellDisStructured(
205 layer=[3, 2, 1],
206 row=[1, 2, 3],
207 column=[2, 2, 2],
208 rate=[-5.0] * 3,
209 ),
210 imod.mf6.MassSourceLoading(
211 rate=get_grid_da(False, np.float64, 0.33),
212 print_input=True,
213 print_flows=False,
214 save_flows=False,
215 ),
216] + [
217 *create_instance_packages(is_unstructured=False),
218 *create_instance_boundary_condition_packages(False),
219]
222UNSTRUCTURED_GRID_PACKAGES = (
223 [
224 imod.mf6.WellDisVertices(
225 layer=[1, 2, 1],
226 cell2d=[3, 12, 23],
227 rate=[-0.1, 0.2, 0.3],
228 print_input=False,
229 print_flows=False,
230 save_flows=False,
231 )
232 ]
233 + [create_vertices_discretization()]
234 + [
235 *create_instance_packages(is_unstructured=True),
236 *create_instance_boundary_condition_packages(True),
237 ]
238)
240GRIDLESS_PACKAGES = [
241 imod.mf6.adv.Advection("upstream"),
242 imod.mf6.Buoyancy(
243 reference_density=1000.0,
244 reference_concentration=[4.0, 25.0],
245 density_concentration_slope=[0.7, -0.375],
246 modelname=["gwt-1", "gwt-2"],
247 species=["salinity", "temperature"],
248 ),
249 imod.mf6.OutputControl(),
250 imod.mf6.SolutionPresetSimple(modelnames=["gwf-1"]),
251 imod.mf6.TimeDiscretization(
252 xr.DataArray(
253 data=[0.001, 7.0, 365.0],
254 coords={"time": pd.date_range("2000-01-01", "2000-01-03")},
255 dims=["time"],
256 ),
257 23,
258 1.02,
259 ),
260 imod.mf6.Well(
261 x=[1.0, 6002.0],
262 y=[3.0, 5004.0],
263 screen_top=[0.0, 0.0],
264 screen_bottom=[-10.0, -10.0],
265 rate=[1.0, 3.0],
266 print_flows=True,
267 validate=True,
268 ),
269 imod.mf6.HorizontalFlowBarrierHydraulicCharacteristic(
270 geometry=gpd.GeoDataFrame(
271 geometry=[
272 shapely.linestrings([-1000.0, 1000.0], [500.0, 500.0]),
273 ],
274 data={
275 "hydraulic_characteristic": [1e-3],
276 "ztop": [10.0],
277 "zbottom": [0.0],
278 },
279 ),
280 ),
281 imod.mf6.HorizontalFlowBarrierMultiplier(
282 geometry=gpd.GeoDataFrame(
283 geometry=[
284 shapely.linestrings([-1000.0, 1000.0], [500.0, 500.0]),
285 ],
286 data={
287 "multiplier": [1.5],
288 "ztop": [10.0],
289 "zbottom": [0.0],
290 },
291 ),
292 ),
293 imod.mf6.HorizontalFlowBarrierResistance(
294 geometry=gpd.GeoDataFrame(
295 geometry=[
296 shapely.linestrings([-1000.0, 1000.0], [500.0, 500.0]),
297 ],
298 data={
299 "resistance": [1e3],
300 "ztop": [10.0],
301 "zbottom": [0.0],
302 },
303 ),
304 ),
305 imod.mf6.GWFGWT("model_name1", "model_name2"),
306]
309def create_exchange_package() -> list[GWFGWF]:
310 cell_id1 = xr.DataArray([[1, 1], [2, 1], [3, 1]])
311 cell_id2 = xr.DataArray([[1, 2], [2, 2], [3, 2]])
312 layer = np.array([12, 13, 14])
313 cl1 = np.ones(len(cell_id1))
314 cl2 = np.ones(len(cell_id1))
315 hwva = cl1 + cl2
317 return [
318 imod.mf6.GWFGWF(
319 "submodel_1",
320 "submodel_2",
321 cell_id1=cell_id1,
322 cell_id2=cell_id2,
323 layer=layer,
324 cl1=cl1,
325 cl2=cl2,
326 hwva=hwva,
327 )
328 ]
331ALL_PACKAGE_INSTANCES = (
332 GRIDLESS_PACKAGES
333 + STRUCTURED_GRID_PACKAGES
334 + UNSTRUCTURED_GRID_PACKAGES
335 + create_exchange_package()
336)