Coverage for C:\src\imod-python\imod\mf6\mf6_hfb_adapter.py: 100%
41 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
1from copy import deepcopy
2from typing import Union
4import numpy as np
5import xarray as xr
7from imod.mf6.boundary_condition import BoundaryCondition
8from imod.mf6.write_context import WriteContext
9from imod.schemata import (
10 AnyValueSchema,
11 CoordsSchema,
12 DimsSchema,
13 DTypeSchema,
14 IndexesSchema,
15)
18class Mf6HorizontalFlowBarrier(BoundaryCondition):
19 """
20 Horizontal Flow Barrier (HFB) package
22 Input to the Horizontal Flow Barrier (HFB) Package is read from the file
23 that has type "HFB6" in the Name File. Only one HFB Package can be
24 specified for a GWF model.
25 https://water.usgs.gov/water-resources/software/MODFLOW-6/mf6io_6.2.2.pdf
27 We recommend using either the ``HorizontalFlowBarrierResistance`` or
28 ``HorizontalFlowBarrierMultiplier`` over this class.
30 Parameters
31 ----------
32 cell_id1: xr.DataArray
33 the cell id on 1 side of the barrier. This is either in column, row format for structured grids or in cell2d
34 format for unstructured grids. The DataArray should contain a coordinate "cell_dims1" which specifies the
35 cell_id structure eg:
36 cell_dims1 (cell_dims1) <U8 'row_1' 'column_1' or cell_dims1 (cell_dims1) <U8 'cell2d_1' 'cell2d_1'
37 cell_id2: xr.DataArray
38 the cell id on the other side of the barrier. This is either in column, row format for structured grids or in
39 cell2d format for unstructured grids. The DataArray should contain a coordinate "cell_dims2" which specifies the
40 cell_id structure eg:
41 cell_dims2 (cell_dims2) <U8 'row_2' 'column_2' or cell_dims2 (cell_dims2) <U8 'cell2d_2' 'cell2d_2'
42 layer: xr.DataArray
43 the layers in which the barrier is active
44 hydraulic_characteristic: xr.DataArray
45 hydraulic characteristic of the barrier: the inverse of the hydraulic
46 resistance. Negative values are interpreted as a multiplier of the cell
47 to cell conductance.
49 Examples
50 --------
51 >> # Structured grid:
52 >> row_1 = [1, 2]
53 >> column_1 = [1, 1]
54 >> row_2 = [1, 2]
55 >> column_2 = [2, 2]
56 >> layer = [1, 2, 3]
57 >>
58 >> cell_indices = np.arange(len(row_1)) + 1
59 >>
60 >> barrier = xr.Dataset()
61 >> barrier["cell_id1"] = xr.DataArray([row_1, column_1], coords={"cell_idx": cell_indices, "cell_dims1": ["row_1", "column_1"]})
62 >> barrier["cell_id2"] = xr.DataArray([row_2, column_2], coords={"cell_idx": cell_indices, "cell_dims2": ["row_2", "column_2"]})
63 >> barrier["hydraulic_characteristic"] = xr.DataArray(np.full((len(layer), len(cell_indices)), 1e-3), coords={"layer": layer, "cell_idx": cell_indices})
64 >> barrier = (barrier.stack(cell_id=("layer", "cell_idx"), create_index=False).drop_vars("cell_idx").reset_coords())
65 >>
66 >> Mf6HorizontalFlowBarrier(**ds)
67 >>
68 >>
69 >> # Unstructured grid
70 >> cell2d_id1 = [1, 2]
71 >> cell2d_id2 = [3, 4]
72 >> layer = [1, 2, 3]
73 >>
74 >> cell_indices = np.arange(len(cell2d_id1)) + 1
75 >>
76 >> barrier = xr.Dataset()
77 >> barrier["cell_id1"] = xr.DataArray(np.array([cell2d_id1]).T, coords={"cell_idx": cell_indices, "cell_dims1": ["cell2d_1"]})
78 >> barrier["cell_id2"] = xr.DataArray(np.array([cell2d_id2]).T, coords={"cell_idx": cell_indices, "cell_dims2": ["cell2d_2"]})
79 >> barrier["hydraulic_characteristic"] = xr.DataArray(np.full((len(layer), len(cell_indices)), 1e-3), coords={"layer": layer, "cell_idx": cell_indices})>>
80 >> barrier = (barrier.stack(cell_id=("layer", "cell_idx"), create_index=False).drop_vars("cell_idx").reset_coords())
81 >>
82 >> Mf6HorizontalFlowBarrier(**ds)
83 >>
85 """
87 _pkg_id = "hfb"
88 _init_schemata = {
89 "hydraulic_characteristic": [
90 DTypeSchema(np.floating),
91 IndexesSchema(),
92 CoordsSchema(("layer",)),
93 DimsSchema("layer", "{edge_dim}") | DimsSchema("{edge_dim}"),
94 ],
95 "idomain": [
96 DTypeSchema(np.integer),
97 DimsSchema("idomain_layer", "y", "x")
98 | DimsSchema("idomain_layer", "{face_dim}"),
99 IndexesSchema(),
100 ],
101 }
102 _write_schemata = {
103 "idomain": (AnyValueSchema(">", 0),),
104 }
105 _period_data = ("hydraulic_characteristic",)
106 _keyword_map = {}
107 _template = BoundaryCondition._initialize_template(_pkg_id)
109 def __init__(
110 self,
111 cell_id1: xr.DataArray,
112 cell_id2: xr.DataArray,
113 layer: xr.DataArray,
114 hydraulic_characteristic: xr.DataArray,
115 print_input: Union[bool, xr.DataArray] = False,
116 validate: Union[bool, xr.DataArray] = True,
117 ):
118 dict_dataset = {
119 "cell_id1": cell_id1,
120 "cell_id2": cell_id2,
121 "layer": layer,
122 "hydraulic_characteristic": hydraulic_characteristic,
123 "print_input": print_input,
124 }
125 super().__init__(dict_dataset)
127 def _get_bin_ds(self):
128 bin_ds = self.dataset[
129 ["cell_id1", "cell_id2", "layer", "hydraulic_characteristic"]
130 ]
132 return bin_ds
134 def _ds_to_arrdict(self, ds):
135 arrdict = super()._ds_to_arrdict(ds)
136 arrdict["cell_dims1"] = ds.coords["cell_dims1"].values
137 arrdict["cell_dims2"] = ds.coords["cell_dims2"].values
139 return arrdict
141 def _to_struct_array(self, arrdict, layer):
142 field_spec = (
143 [("layer_1", np.int32)]
144 + [(dim, np.int32) for dim in arrdict["cell_dims1"]]
145 + [("layer_2", np.int32)]
146 + [(dim, np.int32) for dim in arrdict["cell_dims2"]]
147 + [("hydraulic_characteristic", np.float64)]
148 )
150 sparse_dtype = np.dtype(field_spec)
152 recarr = np.empty(len(arrdict["hydraulic_characteristic"]), dtype=sparse_dtype)
153 recarr["layer_1"] = arrdict["layer"]
154 for dim, value in dict(zip(arrdict["cell_dims1"], arrdict["cell_id1"])).items():
155 recarr[dim] = value
156 recarr["layer_2"] = arrdict["layer"]
157 for dim, value in dict(zip(arrdict["cell_dims2"], arrdict["cell_id2"])).items():
158 recarr[dim] = value
159 recarr["hydraulic_characteristic"] = arrdict["hydraulic_characteristic"]
161 return recarr
163 def write(
164 self,
165 pkgname: str,
166 globaltimes: Union[list[np.datetime64], np.ndarray],
167 write_context: WriteContext,
168 ):
169 # MODFLOW6 does not support binary HFB input.
170 hfb_write_context = deepcopy(write_context)
171 hfb_write_context.use_binary = False
172 super().write(pkgname, globaltimes, hfb_write_context)