Coverage for C:\src\imod-python\imod\mf6\mf6_hfb_adapter.py: 100%

41 statements  

« prev     ^ index     » next       coverage.py v7.4.4, created at 2024-04-08 13:27 +0200

1from copy import deepcopy 

2from typing import Union 

3 

4import numpy as np 

5import xarray as xr 

6 

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) 

16 

17 

18class Mf6HorizontalFlowBarrier(BoundaryCondition): 

19 """ 

20 Horizontal Flow Barrier (HFB) package 

21 

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 

26 

27 We recommend using either the ``HorizontalFlowBarrierResistance`` or 

28 ``HorizontalFlowBarrierMultiplier`` over this class. 

29 

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. 

48 

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 >> 

84 

85 """ 

86 

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) 

108 

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) 

126 

127 def _get_bin_ds(self): 

128 bin_ds = self.dataset[ 

129 ["cell_id1", "cell_id2", "layer", "hydraulic_characteristic"] 

130 ] 

131 

132 return bin_ds 

133 

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 

138 

139 return arrdict 

140 

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 ) 

149 

150 sparse_dtype = np.dtype(field_spec) 

151 

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"] 

160 

161 return recarr 

162 

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)