1import subprocess
2from pathlib import Path
3
4import numpy as np
5import xarray as xr
6
7import imod
8
9nlay = 3
10nrow = 15
11ncol = 15
12shape = (nlay, nrow, ncol)
13
14dx = 5000.0
15dy = -5000.0
16xmin = 0.0
17xmax = dx * ncol
18ymin = 0.0
19ymax = abs(dy) * nrow
20dims = ("layer", "y", "x")
21
22layer = np.array([1, 2, 3])
23y = np.arange(ymax, ymin, dy) + 0.5 * dy
24x = np.arange(xmin, xmax, dx) + 0.5 * dx
25coords = {"layer": layer, "y": y, "x": x}
26
27# Discretization data
28idomain = xr.DataArray(np.ones(shape), coords=coords, dims=dims)
29bottom = xr.DataArray([-200.0, -300.0, -450.0], {"layer": layer}, ("layer",))
30
31# Constant head
32head = xr.full_like(idomain, np.nan).sel(layer=[1, 2])
33head[...] = np.nan
34head[..., 0] = 0.0
35
36# Drainage
37elevation = xr.full_like(idomain.sel(layer=1), np.nan)
38conductance = xr.full_like(idomain.sel(layer=1), np.nan)
39elevation[7, 1:10] = np.array([0.0, 0.0, 10.0, 20.0, 30.0, 50.0, 70.0, 90.0, 100.0])
40conductance[7, 1:10] = 1.0
41
42# Node properties
43icelltype = xr.DataArray([1, 0, 0], {"layer": layer}, ("layer",))
44k = xr.DataArray([1.0e-3, 1.0e-4, 2.0e-4], {"layer": layer}, ("layer",))
45k33 = xr.DataArray([2.0e-8, 2.0e-8, 2.0e-8], {"layer": layer}, ("layer",))
46
47# Recharge
48rch_rate = xr.full_like(idomain.sel(layer=1), 3.0e-8)
49
50# Well
51layer = [3, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
52row = [5, 4, 6, 9, 9, 9, 9, 11, 11, 11, 11, 13, 13, 13, 13]
53column = [11, 6, 12, 8, 10, 12, 14, 8, 10, 12, 14, 8, 10, 12, 14]
54rate = [-5.0] * 15
55
56# Create and fill the groundwater model.
57gwf_model = imod.mf6.GroundwaterFlowModel()
58gwf_model["dis"] = imod.mf6.StructuredDiscretization(
59 top=200.0, bottom=bottom, idomain=idomain
60)
61gwf_model["chd"] = imod.mf6.ConstantHead(
62 head, print_input=True, print_flows=True, save_flows=True
63)
64gwf_model["drn"] = imod.mf6.Drainage(
65 elevation=elevation,
66 conductance=conductance,
67 print_input=True,
68 print_flows=True,
69 save_flows=True,
70)
71gwf_model["ic"] = imod.mf6.InitialConditions(head=0.0)
72gwf_model["npf"] = imod.mf6.NodePropertyFlow(
73 icelltype=icelltype,
74 k=k,
75 k33=k33,
76 variable_vertical_conductance=True,
77 dewatered=True,
78 perched=True,
79 save_flows=True,
80)
81gwf_model["oc"] = imod.mf6.OutputControl(save_head="all", save_budget="all")
82gwf_model["rch"] = imod.mf6.Recharge(rch_rate)
83gwf_model["wel"] = imod.mf6.Well(
84 layer=layer,
85 row=row,
86 column=column,
87 rate=rate,
88 print_input=True,
89 print_flows=True,
90 save_flows=True,
91)
92
93# Attach it to a simulation
94simulation = imod.mf6.Modflow6Simulation("ex01-twri")
95simulation["GWF_1"] = gwf_model
96# Define solver settings
97simulation["solver"] = imod.mf6.Solution(
98 print_option="summary",
99 csv_output=False,
100 no_ptc=True,
101 outer_hclose=1.0e-4,
102 outer_maximum=500,
103 under_relaxation=None,
104 inner_hclose=1.0e-4,
105 inner_rclose=0.001,
106 inner_maximum=100,
107 linear_acceleration="cg",
108 scaling_method=None,
109 reordering_method=None,
110 relaxation_factor=0.97,
111)
112# Collect time discretization
113simulation.time_discretization(times=["2000-01-01", "2000-01-02"])
114
115modeldir = Path("ex01-twri")
116simulation.write(modeldir)
117with imod.util.cd(modeldir):
118 p = subprocess.run("mf6", check=True, capture_output=True, text=True)
119 assert p.stdout.endswith("Normal termination of simulation.\n")
120
121head = imod.mf6.open_hds(modeldir / "GWF_1/GWF_1.hds", modeldir / "GWF_1/dis.dis.grb")