Coverage for fixtures\mf6_rectangle_with_lakes.py: 100%

56 statements  

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

1import numpy as np 

2import pytest 

3import xarray as xr 

4 

5import imod 

6from imod.mf6.lak import Lake, LakeData 

7 

8 

9def create_lake(idomain, xmin_index, xmax_index, ymin_index, ymax_index, name): 

10 is_lake = xr.full_like(idomain, False, dtype=bool) 

11 is_lake.values[0, xmin_index : xmax_index + 1, ymin_index : ymax_index + 1] = True 

12 

13 lake_table = None 

14 return create_lake_data_structured( 

15 is_lake, starting_stage=11.0, name=name, lake_table=lake_table 

16 ) 

17 

18 

19@pytest.fixture(scope="function") 

20def rectangle_with_lakes(): 

21 shape = nlay, nrow, ncol = 2, 30, 30 

22 

23 dx = 100.0 

24 dy = -100.0 

25 

26 xmin = 0.0 

27 xmax = dx * ncol 

28 ymin = 0.0 

29 ymax = abs(dy) * nrow 

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

31 layer = np.arange(1, nlay + 1) 

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

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

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

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

36 

37 # Discretization data 

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

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

40 idomain = like.astype(np.int8) 

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

42 k33 = xr.DataArray([2.0e-4, 1.0e-4], {"layer": layer}, ("layer",)) 

43 

44 # Constant head 

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

46 head[..., 0] = 0.0 

47 

48 # Create and fill the groundwater model. 

49 gwf_model = imod.mf6.GroundwaterFlowModel() 

50 gwf_model["dis"] = imod.mf6.StructuredDiscretization( 

51 top=0.0, bottom=bottom, idomain=idomain 

52 ) 

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

54 head, print_input=True, print_flows=True, save_flows=True 

55 ) 

56 gwf_model["chd"] = imod.mf6.SpecificStorage( 

57 specific_storage=1.0e-5, 

58 specific_yield=0.15, 

59 convertible=0, 

60 transient=True, 

61 ) 

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

63 

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

65 icelltype=icelltype, 

66 k=k, 

67 k33=k33, 

68 variable_vertical_conductance=True, 

69 dewatered=True, 

70 perched=True, 

71 save_flows=True, 

72 ) 

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

74 

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

76 simulation["GWF_1"] = gwf_model 

77 # Define solver settings 

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

79 modelnames=["GWF_1"], 

80 print_option="summary", 

81 csv_output=False, 

82 no_ptc=True, 

83 outer_dvclose=1.0e-4, 

84 outer_maximum=500, 

85 under_relaxation=None, 

86 inner_dvclose=1.0e-4, 

87 inner_rclose=0.001, 

88 inner_maximum=100, 

89 linear_acceleration="cg", 

90 scaling_method=None, 

91 reordering_method=None, 

92 relaxation_factor=0.97, 

93 ) 

94 # Collect time discretization 

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

96 lake1 = create_lake(idomain, 2, 3, 2, 3, "first_lake") 

97 lake2 = create_lake(idomain, 12, 13, 12, 13, "second_lake") 

98 simulation["GWF_1"]["lake"] = Lake.from_lakes_and_outlets([lake1, lake2]) 

99 return simulation 

100 

101 

102def create_lake_data_structured( 

103 is_lake, 

104 starting_stage, 

105 name, 

106 status=None, 

107 stage=None, 

108 rainfall=None, 

109 evaporation=None, 

110 runoff=None, 

111 inflow=None, 

112 withdrawal=None, 

113 auxiliary=None, 

114 lake_table=None, 

115): 

116 HORIZONTAL = 0 

117 connection_type = xr.full_like(is_lake, HORIZONTAL, dtype=np.float64).where(is_lake) 

118 bed_leak = xr.full_like(is_lake, 0.2, dtype=np.float64).where(is_lake) 

119 top_elevation = xr.full_like(is_lake, 0.0, dtype=np.float64).where(is_lake) 

120 bot_elevation = xr.full_like(is_lake, -1.0, dtype=np.float64).where(is_lake) 

121 connection_length = xr.full_like(is_lake, 0.5, dtype=np.float64).where(is_lake) 

122 connection_width = xr.full_like(is_lake, 0.6, dtype=np.float64).where(is_lake) 

123 return LakeData( 

124 starting_stage=starting_stage, 

125 boundname=name, 

126 connection_type=connection_type, 

127 bed_leak=bed_leak, 

128 top_elevation=top_elevation, 

129 bot_elevation=bot_elevation, 

130 connection_length=connection_length, 

131 connection_width=connection_width, 

132 status=status, 

133 stage=stage, 

134 rainfall=rainfall, 

135 evaporation=evaporation, 

136 runoff=runoff, 

137 inflow=inflow, 

138 withdrawal=withdrawal, 

139 auxiliary=auxiliary, 

140 lake_table=lake_table, 

141 )