Coverage for fixtures\mf6_twri_disv_fixture.py: 100%

65 statements  

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

1import numpy as np 

2import pytest 

3import xarray as xr 

4import xugrid as xu 

5 

6import imod 

7 

8 

9def make_twri_disv_model(): 

10 nlay = 3 

11 nrow = 15 

12 ncol = 15 

13 shape = (nlay, nrow, ncol) 

14 

15 dx = 5000.0 

16 dy = -5000.0 

17 xmin = 0.0 

18 xmax = dx * ncol 

19 ymin = 0.0 

20 ymax = abs(dy) * nrow 

21 dims = ("layer", "y", "x") 

22 

23 layer = np.array([1, 2, 3]) 

24 y = np.arange(ymax, ymin, dy) + 0.5 * dy 

25 x = np.arange(xmin, xmax, dx) + 0.5 * dx 

26 coords = {"layer": layer, "y": y, "x": x} 

27 

28 # Discretization data 

29 like = xr.DataArray(np.ones(shape), coords=coords, dims=dims) 

30 idomain = like.astype(np.int8) 

31 

32 bottom = xr.DataArray([-200.0, -300.0, -450.0], {"layer": layer}, ("layer",)) 

33 

34 # Constant head 

35 head = xr.full_like(like, np.nan).sel(layer=[1, 2]) 

36 head[..., 0] = 0.0 

37 

38 # Drainage 

39 elevation = xr.full_like(like.sel(layer=1), np.nan) 

40 conductance = xr.full_like(like.sel(layer=1), np.nan) 

41 elevation[7, 1:10] = np.array([0.0, 0.0, 10.0, 20.0, 30.0, 50.0, 70.0, 90.0, 100.0]) 

42 conductance[7, 1:10] = 1.0 

43 

44 # Node properties 

45 icelltype = xr.DataArray([1, 0, 0], {"layer": layer}, ("layer",)) 

46 k = xr.DataArray([1.0e-3, 1.0e-4, 2.0e-4], {"layer": layer}, ("layer",)) 

47 k33 = xr.DataArray([2.0e-8, 2.0e-8, 2.0e-8], {"layer": layer}, ("layer",)) 

48 

49 # Recharge 

50 rch_rate = xr.full_like(like.sel(layer=1), 3.0e-8) 

51 

52 # Well 

53 layer = np.array([3, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]) 

54 row = np.array([5, 4, 6, 9, 9, 9, 9, 11, 11, 11, 11, 13, 13, 13, 13]) 

55 column = np.array([11, 6, 12, 8, 10, 12, 14, 8, 10, 12, 14, 8, 10, 12, 14]) 

56 rate = np.array( 

57 [ 

58 -5.0, 

59 -5.0, 

60 -5.0, 

61 -5.0, 

62 -5.0, 

63 -5.0, 

64 -5.0, 

65 -5.0, 

66 -5.0, 

67 -5.0, 

68 -5.0, 

69 -5.0, 

70 -5.0, 

71 -5.0, 

72 -5.0, 

73 ] 

74 ) 

75 cell2d = (row - 1) * ncol + column 

76 

77 # %% 

78 

79 grid = xu.Ugrid2d.from_structured(idomain) 

80 

81 def toface(da): 

82 return imod.util.spatial.ugrid2d_data(da, grid.face_dimension) 

83 

84 idomain = xu.UgridDataArray(toface(idomain), grid) 

85 head = xu.UgridDataArray(toface(head), grid) 

86 elevation = xu.UgridDataArray(toface(elevation), grid).assign_coords(layer=1) 

87 conductance = xu.UgridDataArray(toface(conductance), grid).assign_coords(layer=1) 

88 rch_rate = xu.UgridDataArray(toface(rch_rate), grid).assign_coords(layer=1) 

89 

90 # Create and fill the groundwater model. 

91 gwf_model = imod.mf6.GroundwaterFlowModel() 

92 gwf_model["dis"] = imod.mf6.VerticesDiscretization( 

93 top=200.0, bottom=bottom, idomain=idomain 

94 ) 

95 gwf_model["chd"] = imod.mf6.ConstantHead( 

96 head, print_input=True, print_flows=True, save_flows=True 

97 ) 

98 gwf_model["drn"] = imod.mf6.Drainage( 

99 elevation=elevation, 

100 conductance=conductance, 

101 print_input=True, 

102 print_flows=True, 

103 save_flows=True, 

104 ) 

105 gwf_model["ic"] = imod.mf6.InitialConditions(start=0.0) 

106 gwf_model["npf"] = imod.mf6.NodePropertyFlow( 

107 icelltype=icelltype, 

108 k=k, 

109 k33=k33, 

110 variable_vertical_conductance=True, 

111 dewatered=True, 

112 perched=True, 

113 save_flows=True, 

114 ) 

115 gwf_model["oc"] = imod.mf6.OutputControl(save_head="all", save_budget="all") 

116 gwf_model["rch"] = imod.mf6.Recharge(rch_rate) 

117 gwf_model["sto"] = imod.mf6.SpecificStorage( 

118 specific_storage=1.0e-15, 

119 specific_yield=0.15, 

120 convertible=0, 

121 transient=False, 

122 ) 

123 gwf_model["wel"] = imod.mf6.WellDisVertices( 

124 layer=layer, 

125 cell2d=cell2d, 

126 rate=rate, 

127 print_input=True, 

128 print_flows=True, 

129 save_flows=True, 

130 ) 

131 

132 # Attach it to a simulation 

133 simulation = imod.mf6.Modflow6Simulation("ex01-twri-disv") 

134 simulation["GWF_1"] = gwf_model 

135 # Define solver settings 

136 simulation["solver"] = imod.mf6.Solution( 

137 modelnames=["GWF_1"], 

138 print_option="summary", 

139 outer_dvclose=1.0e-4, 

140 outer_maximum=500, 

141 under_relaxation=None, 

142 inner_dvclose=1.0e-4, 

143 inner_rclose=0.001, 

144 inner_maximum=100, 

145 linear_acceleration="cg", 

146 scaling_method=None, 

147 reordering_method=None, 

148 relaxation_factor=0.97, 

149 ) 

150 # Collect time discretization 

151 simulation.create_time_discretization(additional_times=["2000-01-01", "2000-01-02"]) 

152 

153 return simulation 

154 

155 

156@pytest.fixture(scope="session") 

157def twri_disv_model(): 

158 return make_twri_disv_model()