Coverage for fixtures\flow_example_fixture.py: 11%
70 statements
« prev ^ index » next coverage.py v7.5.1, created at 2024-05-08 14:15 +0200
« prev ^ index » next coverage.py v7.5.1, created at 2024-05-08 14:15 +0200
1import numpy as np
2import pandas as pd
3import pytest
4import xarray as xr
6import imod.flow as flow
9def create_imodflow_model():
10 shape = nlay, nrow, ncol = 3, 9, 9
12 dx = 100.0
13 dy = -100.0
14 dz = np.array([5, 30, 100])
15 xmin = 0.0
16 xmax = dx * ncol
17 ymin = 0.0
18 ymax = abs(dy) * nrow
19 dims = ("layer", "y", "x")
20 kh = 10.0
21 kva = 1.0
22 sto = 0.001
24 # Some coordinates
25 times = pd.date_range(start="1/1/2018", end="12/1/2018", freq="MS")
26 layer = np.arange(1, nlay + 1)
27 y = np.arange(ymax, ymin, dy) + 0.5 * dy
28 x = np.arange(xmin, xmax, dx) + 0.5 * dx
29 coords = {"layer": layer, "y": y, "x": x}
31 # top & bottom
32 surface = 0.0
33 interfaces = np.insert((surface - np.cumsum(dz)), 0, surface)
34 bottom = xr.DataArray(interfaces[1:], coords={"layer": layer}, dims="layer")
35 top = xr.DataArray(interfaces[:-1], coords={"layer": layer}, dims="layer")
37 # constant head & ibound
38 ibound = xr.DataArray(np.ones(shape), coords=coords, dims=dims)
39 starting_head = ibound.copy()
40 trend = np.ones(times[:-1].shape)
41 trend = np.cumsum(trend)
42 trend_da = xr.DataArray(trend, coords={"time": times[:-1]}, dims=["time"])
43 is_x_edge = starting_head.x.isin([x[0], x[-1]])
44 head_edge = starting_head.where(is_x_edge)
45 head_edge_rising = trend_da * head_edge
46 ibound = ibound.where(head_edge.isnull(), other=-2)
48 # general head boundary
49 is_x_central = starting_head.x == x[4]
50 head_central = starting_head.where(is_x_central).sel(layer=1)
51 period_times = times[[0, 6]] - np.timedelta64(365, "D")
52 periods_da = xr.DataArray([10, 4], coords={"time": period_times}, dims=["time"])
53 head_periodic = periods_da * head_central
54 timemap = {
55 period_times[0]: "winter",
56 period_times[1]: "summer",
57 }
59 # Wells
60 wel_df = pd.DataFrame()
61 wel_df["id_name"] = np.arange(len(x)).astype(str)
62 wel_df["x"] = x
63 wel_df["y"] = y
64 wel_df["rate"] = dx * dy * -1 * 0.5
65 wel_df["time"] = np.tile(times[:-1], 2)[: len(x)]
66 wel_df["layer"] = 2
68 # Fill the model
69 m = flow.ImodflowModel("imodflow")
70 m["pcg"] = flow.PreconditionedConjugateGradientSolver()
71 m["bnd"] = flow.Boundary(ibound)
72 m["top"] = flow.Top(top)
73 m["bottom"] = flow.Bottom(bottom)
74 m["khv"] = flow.HorizontalHydraulicConductivity(kh)
75 m["kva"] = flow.VerticalAnisotropy(kva)
76 m["sto"] = flow.StorageCoefficient(sto)
77 m["shd"] = flow.StartingHead(starting_head)
78 m["chd"] = flow.ConstantHead(head=head_edge_rising)
79 m["ghb"] = flow.GeneralHeadBoundary(head=head_periodic, conductance=10.0)
80 m["ghb"].periodic_stress(timemap)
81 m["ghb2"] = flow.GeneralHeadBoundary(head=head_periodic + 10.0, conductance=1.0)
82 m["ghb2"].periodic_stress(timemap)
83 m["wel"] = flow.Well(**wel_df)
84 m["oc"] = flow.OutputControl(save_head=-1, save_flux=-1)
85 m.create_time_discretization(additional_times=times[-1])
86 return m
89@pytest.fixture(scope="session")
90def imodflow_model():
91 return create_imodflow_model()