Coverage for C:\src\imod-python\imod\prepare\wells.py: 100%

71 statements  

« prev     ^ index     » next       coverage.py v7.5.1, created at 2024-05-08 14:15 +0200

1""" 

2Assign wells to layers. 

3""" 

4 

5from typing import Optional, Union 

6 

7import numpy as np 

8import pandas as pd 

9import xarray as xr 

10import xugrid as xu 

11 

12import imod 

13 

14 

15def vectorized_overlap(bounds_a, bounds_b): 

16 """ 

17 Vectorized overlap computation. 

18 Compare with: 

19 overlap = max(0, min(a[1], b[1]) - max(a[0], b[0])) 

20 """ 

21 return np.maximum( 

22 0.0, 

23 np.minimum(bounds_a[:, 1], bounds_b[:, 1]) 

24 - np.maximum(bounds_a[:, 0], bounds_b[:, 0]), 

25 ) 

26 

27 

28def compute_overlap(wells, top, bottom): 

29 # layer bounds shape of (n_well, n_layer, 2) 

30 layer_bounds = np.stack((bottom, top), axis=-1) 

31 well_bounds = np.broadcast_to( 

32 np.stack( 

33 (wells["bottom"].to_numpy(), wells["top"].to_numpy()), 

34 axis=-1, 

35 )[np.newaxis, :, :], 

36 layer_bounds.shape, 

37 ) 

38 overlap = vectorized_overlap( 

39 well_bounds.reshape((-1, 2)), 

40 layer_bounds.reshape((-1, 2)), 

41 ) 

42 return overlap 

43 

44 

45def locate_wells( 

46 wells: pd.DataFrame, 

47 top: Union[xr.DataArray, xu.UgridDataArray], 

48 bottom: Union[xr.DataArray, xu.UgridDataArray], 

49 k: Optional[Union[xr.DataArray, xu.UgridDataArray]], 

50 validate: bool = True, 

51): 

52 if not isinstance(top, (xu.UgridDataArray, xr.DataArray)): 

53 raise TypeError( 

54 "top and bottom should be DataArray or UgridDataArray, received: " 

55 f"{type(top).__name__}" 

56 ) 

57 

58 # Default to a xy_k value of 1.0: weigh every layer equally. 

59 xy_k = 1.0 

60 first = wells.groupby("id").first() 

61 x = first["x"].to_numpy() 

62 y = first["y"].to_numpy() 

63 

64 xy_top = imod.select.points_values(top, x=x, y=y, out_of_bounds="ignore") 

65 xy_bottom = imod.select.points_values(bottom, x=x, y=y, out_of_bounds="ignore") 

66 

67 # Raise exception if not all wells could be mapped onto the domain 

68 if validate and len(x) > len(xy_top["index"]): 

69 inside = imod.select.points_in_bounds(top, x=x, y=y) 

70 out = np.where(~inside) 

71 raise ValueError( 

72 f"well at x = {x[out[0]]} and y = {y[out[0]]} could not be mapped on the grid" 

73 ) 

74 

75 if k is not None: 

76 xy_k = imod.select.points_values(k, x=x, y=y, out_of_bounds="ignore") 

77 

78 # Discard out-of-bounds wells. 

79 index = xy_top["index"] 

80 if validate and not np.array_equal(xy_bottom["index"], index): 

81 raise ValueError("bottom grid does not match top grid") 

82 if validate and k is not None and not np.array_equal(xy_k["index"], index): # type: ignore 

83 raise ValueError("k grid does not match top grid") 

84 id_in_bounds = first.index[index] 

85 

86 return id_in_bounds, xy_top, xy_bottom, xy_k 

87 

88 

89def assign_wells( 

90 wells: pd.DataFrame, 

91 top: Union[xr.DataArray, xu.UgridDataArray], 

92 bottom: Union[xr.DataArray, xu.UgridDataArray], 

93 k: Optional[Union[xr.DataArray, xu.UgridDataArray]] = None, 

94 minimum_thickness: Optional[float] = 0.05, 

95 minimum_k: Optional[float] = 1.0, 

96 validate: bool = True, 

97) -> pd.DataFrame: 

98 """ 

99 Distribute well pumping rate according to filter length when ``k=None``, or 

100 to transmissivity of the sediments surrounding the filter. Minimum 

101 thickness and minimum k should be set to avoid placing wells in clay 

102 layers. 

103 

104 Wells located outside of the grid are removed. 

105 

106 Parameters 

107 ---------- 

108 wells: pd.DataFrame 

109 Should contain columns x, y, id, top, bottom, rate. 

110 top: xr.DataArray or xu.UgridDataArray 

111 Top of the model layers. 

112 bottom: xr.DataArray or xu.UgridDataArray 

113 Bottom of the model layers. 

114 k: xr.DataArray or xu.UgridDataArray, optional 

115 Horizontal conductivity of the model layers. 

116 minimum_thickness: float, optional, default: 0.01 

117 minimum_k: float, optional, default: 1.0 

118 Minimum conductivity 

119 validate: bool 

120 raise an excpetion if one of the wells is not in the domain 

121 Returns 

122 ------- 

123 placed_wells: pd.DataFrame 

124 Wells with rate subdivided per layer. Contains the original columns of 

125 ``wells``, as well as layer, overlap, transmissivity. 

126 """ 

127 

128 names = {"x", "y", "id", "top", "bottom", "rate"} 

129 missing = names.difference(wells.columns) 

130 if missing: 

131 raise ValueError(f"Columns are missing in wells dataframe: {missing}") 

132 

133 types = [type(arg) for arg in (top, bottom, k) if arg is not None] 

134 if len(set(types)) != 1: 

135 members = ",".join([t.__name__ for t in types]) 

136 raise TypeError( 

137 "top, bottom, and optionally k should be of the same type, " 

138 f"received: {members}" 

139 ) 

140 

141 id_in_bounds, xy_top, xy_bottom, xy_k = locate_wells( 

142 wells, top, bottom, k, validate 

143 ) 

144 wells_in_bounds = wells.set_index("id").loc[id_in_bounds].reset_index() 

145 first = wells_in_bounds.groupby("id").first() 

146 overlap = compute_overlap(first, xy_top, xy_bottom) 

147 

148 if k is None: 

149 k = 1.0 

150 else: 

151 k = xy_k.values.ravel() 

152 

153 # Distribute rate according to transmissivity. 

154 n_layer, n_well = xy_top.shape 

155 df = pd.DataFrame( 

156 index=pd.Index(np.tile(first.index, n_layer), name="id"), 

157 data={ 

158 "layer": np.repeat(top["layer"], n_well), 

159 "overlap": overlap, 

160 "k": k, 

161 "transmissivity": overlap * k, 

162 }, 

163 ) 

164 # remove entries 

165 # -in very thin layers or when the wellbore penetrates the layer very little 

166 # -in low conductivity layers 

167 df = df.loc[(df["overlap"] >= minimum_thickness) & (df["k"] >= minimum_k)] 

168 df["rate"] = df["transmissivity"] / df.groupby("id")["transmissivity"].transform( 

169 "sum" 

170 ) 

171 # Create a unique index for every id-layer combination. 

172 df["index"] = np.arange(len(df)) 

173 df = df.reset_index() 

174 

175 # Get rid of those that are removed because of minimum thickness or 

176 # transmissivity. 

177 wells_in_bounds = wells_in_bounds.loc[wells_in_bounds["id"].isin(df["id"].unique())] 

178 

179 # Use pandas multi-index broadcasting. 

180 # Maintain all other columns as-is. 

181 wells_in_bounds["index"] = 1 # N.B. integer! 

182 wells_in_bounds["overlap"] = 1.0 

183 wells_in_bounds["k"] = 1.0 

184 wells_in_bounds["transmissivity"] = 1.0 

185 columns = list(set(wells_in_bounds.columns).difference(df.columns)) 

186 

187 indexes = ["id"] 

188 for dim in ["species", "time"]: 

189 if dim in wells_in_bounds: 

190 indexes.append(dim) 

191 columns.remove(dim) 

192 

193 df[columns] = 1 # N.B. integer! 

194 

195 assigned = ( 

196 wells_in_bounds.set_index(indexes) * df.set_index(["id", "layer"]) 

197 ).reset_index() 

198 return assigned