| 1 | """ Testing CULVERT (Changing from Horizontal Abstraction to Vertical Abstraction |
|---|
| 2 | |
|---|
| 3 | This example includes a Model Topography that shows a TYPICAL Headwall Configuration |
|---|
| 4 | |
|---|
| 5 | The aim is to change the Culvert Routine to Model more precisely the abstraction |
|---|
| 6 | from a vertical face. |
|---|
| 7 | |
|---|
| 8 | The inflow must include the impact of Approach velocity. |
|---|
| 9 | Similarly the Outflow has MOMENTUM Not just Up welling as in the Horizontal Style |
|---|
| 10 | abstraction |
|---|
| 11 | |
|---|
| 12 | """ |
|---|
| 13 | print 'Starting.... Importing Modules...' |
|---|
| 14 | #------------------------------------------------------------------------------ |
|---|
| 15 | # Import necessary modules |
|---|
| 16 | #------------------------------------------------------------------------------ |
|---|
| 17 | import anuga |
|---|
| 18 | |
|---|
| 19 | from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross |
|---|
| 20 | from anuga.shallow_water.shallow_water_domain import Domain |
|---|
| 21 | from anuga.shallow_water.forcing import Rainfall, Inflow |
|---|
| 22 | #from anuga.shallow_water.forcing import Reflective_boundary |
|---|
| 23 | #from anuga.shallow_water.forcing import Dirichlet_boundary |
|---|
| 24 | #from anuga.shallow_water.forcing import Transmissive_boundary, Time_boundary |
|---|
| 25 | |
|---|
| 26 | #from anuga.culvert_flows.culvert_class import Culvert_flow |
|---|
| 27 | from anuga.structures.boyd_pipe_operator import Boyd_pipe_operator |
|---|
| 28 | from anuga.structures.boyd_box_operator import Boyd_box_operator |
|---|
| 29 | #from anuga.culvert_flows.culvert_routines import weir_orifice_channel_culvert_model |
|---|
| 30 | from math import pi,pow,sqrt |
|---|
| 31 | |
|---|
| 32 | import numpy as num |
|---|
| 33 | |
|---|
| 34 | |
|---|
| 35 | |
|---|
| 36 | |
|---|
| 37 | #------------------------------------------------------------------------------ |
|---|
| 38 | # Setup computational domain |
|---|
| 39 | #------------------------------------------------------------------------------ |
|---|
| 40 | print 'Setting up domain' |
|---|
| 41 | |
|---|
| 42 | length = 120. #x-Dir |
|---|
| 43 | width = 200. #y-dir |
|---|
| 44 | |
|---|
| 45 | dx = dy = 2.0 # Resolution: Length of subdivisions on both axes |
|---|
| 46 | #dx = dy = .5 # Resolution: Length of subdivisions on both axes |
|---|
| 47 | #dx = dy = .5 # Resolution: Length of subdivisions on both axes |
|---|
| 48 | #dx = dy = .1 # Resolution: Length of subdivisions on both axes |
|---|
| 49 | |
|---|
| 50 | points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy), |
|---|
| 51 | len1=length, len2=width) |
|---|
| 52 | domain = Domain(points, vertices, boundary) |
|---|
| 53 | domain.set_name('Test_Outlet_Ctrl') # Output name |
|---|
| 54 | domain.set_default_order(2) |
|---|
| 55 | domain.H0 = 0.01 |
|---|
| 56 | domain.tight_slope_limiters = 1 |
|---|
| 57 | |
|---|
| 58 | print 'Size', len(domain) |
|---|
| 59 | |
|---|
| 60 | #------------------------------------------------------------------------------ |
|---|
| 61 | # Setup initial conditions |
|---|
| 62 | #------------------------------------------------------------------------------ |
|---|
| 63 | |
|---|
| 64 | def topography(x, y): |
|---|
| 65 | """Set up a weir |
|---|
| 66 | |
|---|
| 67 | A culvert will connect either side |
|---|
| 68 | """ |
|---|
| 69 | # General Slope of Topography |
|---|
| 70 | z=10.0-x/100.0 # % Longitudinal Slope |
|---|
| 71 | |
|---|
| 72 | # NOW Add bits and Pieces to topography |
|---|
| 73 | bank_hgt=10.0 |
|---|
| 74 | bridge_width = 50.0 |
|---|
| 75 | bank_width = 10.0 |
|---|
| 76 | |
|---|
| 77 | us_apron_skew = 1.0 # 1.0 = 1 Length: 1 Width, 2.0 = 2 Length : 1 Width |
|---|
| 78 | us_start_x = 10.0 |
|---|
| 79 | top_start_y = 50.0 |
|---|
| 80 | us_slope = 3.0 #Horiz : Vertic |
|---|
| 81 | ds_slope = 3.0 |
|---|
| 82 | ds_apron_skew = 1.0 # 1.0 = 1 Length: 1 Width, 2.0 = 2 Length : 1 Width |
|---|
| 83 | centre_line_y= top_start_y+bridge_width/2.0 |
|---|
| 84 | |
|---|
| 85 | # CALCULATE PARAMETERS TO FORM THE EMBANKMENT |
|---|
| 86 | us_slope_length = bank_hgt*us_slope |
|---|
| 87 | us_end_x =us_start_x + us_slope_length |
|---|
| 88 | us_toe_start_y =top_start_y - us_slope_length / us_apron_skew |
|---|
| 89 | us_toe_end_y = top_start_y + bridge_width + us_slope_length / us_apron_skew |
|---|
| 90 | |
|---|
| 91 | top_end_y = top_start_y + bridge_width |
|---|
| 92 | ds_slope_length = bank_hgt*ds_slope |
|---|
| 93 | ds_start_x = us_end_x + bank_width |
|---|
| 94 | ds_end_x = ds_start_x + ds_slope_length |
|---|
| 95 | |
|---|
| 96 | ds_toe_start_y =top_start_y - ds_slope_length / ds_apron_skew |
|---|
| 97 | ds_toe_end_y = top_start_y + bridge_width + ds_slope_length / ds_apron_skew |
|---|
| 98 | |
|---|
| 99 | |
|---|
| 100 | N = len(x) |
|---|
| 101 | for i in range(N): |
|---|
| 102 | |
|---|
| 103 | # Sloping Embankment Across Channel |
|---|
| 104 | if us_start_x < x[i] < us_end_x +0.1: # For UPSLOPE on the Upstream FACE |
|---|
| 105 | #if 5.0 < x[i] < 10.1: # For a Range of X, and over a Range of Y based on X adjust Z |
|---|
| 106 | if us_toe_start_y +(x[i] - us_start_x)/us_apron_skew < y[i] < us_toe_end_y - (x[i] - us_start_x)/us_apron_skew: |
|---|
| 107 | #if 49.0+(x[i]-5.0)/5.0 < y[i] < 151.0 - (x[i]-5.0)/5.0: # Cut Out Base Segment for Culvert FACE |
|---|
| 108 | z[i]=z[i] # Flat Apron |
|---|
| 109 | #z[i] += z[i] + (x[i] - us_start_x)/us_slope |
|---|
| 110 | #pass |
|---|
| 111 | else: |
|---|
| 112 | z[i] += z[i] + (x[i] - us_start_x)/us_slope # Sloping Segment U/S Face |
|---|
| 113 | if us_end_x < x[i] < ds_start_x + 0.1: |
|---|
| 114 | z[i] += z[i]+bank_hgt # Flat Crest of Embankment |
|---|
| 115 | if ds_start_x < x[i] < ds_end_x: # DOWN SDLOPE Segment on Downstream face |
|---|
| 116 | if top_start_y-(x[i]-ds_start_x)/ds_apron_skew < y[i] < top_end_y + (x[i]-ds_start_x)/ds_apron_skew: # Cut Out Segment for Culvert FACE |
|---|
| 117 | z[i]=z[i] # Flat Apron |
|---|
| 118 | #z[i] += z[i]+bank_hgt-(x[i] -ds_start_x)/ds_slope |
|---|
| 119 | #pass |
|---|
| 120 | else: |
|---|
| 121 | z[i] += z[i]+bank_hgt-(x[i] -ds_start_x)/ds_slope # Sloping D/S Face |
|---|
| 122 | |
|---|
| 123 | |
|---|
| 124 | |
|---|
| 125 | return z |
|---|
| 126 | |
|---|
| 127 | print 'Setting Quantities....' |
|---|
| 128 | domain.set_quantity('elevation', topography) # Use function for elevation |
|---|
| 129 | domain.set_quantity('friction', 0.01) # Constant friction |
|---|
| 130 | domain.set_quantity('stage', |
|---|
| 131 | expression='elevation') # Dry initial condition |
|---|
| 132 | |
|---|
| 133 | |
|---|
| 134 | |
|---|
| 135 | |
|---|
| 136 | #------------------------------------------------------------------------------ |
|---|
| 137 | # Setup specialised forcing terms |
|---|
| 138 | #------------------------------------------------------------------------------ |
|---|
| 139 | |
|---|
| 140 | #------------------------------------------------------------------------------ |
|---|
| 141 | # Setup CULVERT INLETS and OUTLETS in Current Topography |
|---|
| 142 | #------------------------------------------------------------------------------ |
|---|
| 143 | print 'Defining Structures' |
|---|
| 144 | |
|---|
| 145 | # DEFINE CULVERT INLET AND OUTLETS |
|---|
| 146 | |
|---|
| 147 | |
|---|
| 148 | #culvert0 = Culvert_operator(domain, |
|---|
| 149 | # end_point0=[40.0, 75.0], |
|---|
| 150 | # end_point1=[50.0, 75.0], |
|---|
| 151 | # width=50.0, |
|---|
| 152 | # depth=10.0, |
|---|
| 153 | # apron=5.0, |
|---|
| 154 | # verbose=False) |
|---|
| 155 | |
|---|
| 156 | |
|---|
| 157 | #------------------------------------------------------------------------------ |
|---|
| 158 | # Setup culverts |
|---|
| 159 | #------------------------------------------------------------------------------ |
|---|
| 160 | |
|---|
| 161 | culverts = [] |
|---|
| 162 | number_of_culverts = 2 |
|---|
| 163 | for i in range(number_of_culverts): |
|---|
| 164 | culvert_width = 50.0/number_of_culverts |
|---|
| 165 | y = 100-i*culvert_width - culvert_width/2.0 |
|---|
| 166 | ep0 = num.array([40.0, y]) |
|---|
| 167 | ep1 = num.array([50.0, y]) |
|---|
| 168 | |
|---|
| 169 | losses = {'inlet':0.5, 'outlet':1, 'bend':0, 'grate':0, 'pier': 0, 'other': 0} |
|---|
| 170 | # culverts.append(Boyd_pipe_operator(domain, |
|---|
| 171 | # end_point0=ep0, |
|---|
| 172 | # end_point1=ep1, |
|---|
| 173 | # losses=losses, |
|---|
| 174 | # diameter=1.5, #culvert_width, #3.658, |
|---|
| 175 | # apron=6.0, |
|---|
| 176 | # use_momentum_jet=True, |
|---|
| 177 | # use_velocity_head=True, |
|---|
| 178 | # manning=0.013, |
|---|
| 179 | # logging=True, |
|---|
| 180 | # label='culvert', |
|---|
| 181 | # verbose=False)) |
|---|
| 182 | |
|---|
| 183 | |
|---|
| 184 | |
|---|
| 185 | |
|---|
| 186 | |
|---|
| 187 | Boyd_box_operator(domain, |
|---|
| 188 | losses=losses, |
|---|
| 189 | width=culvert_width, |
|---|
| 190 | height=10.0, |
|---|
| 191 | end_points=[ep0,ep1], |
|---|
| 192 | apron=6.0, |
|---|
| 193 | manning=0.013, |
|---|
| 194 | label='box_culvert', |
|---|
| 195 | logging=True, |
|---|
| 196 | verbose=False) |
|---|
| 197 | |
|---|
| 198 | #losses = {'inlet':1, 'outlet':1, 'bend':1, 'grate':1, 'pier': 1, 'other': 1} |
|---|
| 199 | #culvert2 = Culvert_operator(domain, |
|---|
| 200 | #end_point0=[40.0, 62.5], |
|---|
| 201 | #end_point1=[50.0, 62.5], |
|---|
| 202 | #losses, |
|---|
| 203 | #width=25.0, |
|---|
| 204 | #depth=10.0, |
|---|
| 205 | #apron=5.0, |
|---|
| 206 | #manning=0.013, |
|---|
| 207 | #verbose=False) |
|---|
| 208 | |
|---|
| 209 | |
|---|
| 210 | |
|---|
| 211 | #------------------------------------------------------------------------------ |
|---|
| 212 | # Setup boundary conditions |
|---|
| 213 | #------------------------------------------------------------------------------ |
|---|
| 214 | print 'Setting Boundary Conditions' |
|---|
| 215 | Br = anuga.Reflective_boundary(domain) # Solid reflective wall |
|---|
| 216 | Bi = anuga.Dirichlet_boundary([0.0, 0.0, 0.0]) # Inflow based on Flow Depth and Approaching Momentum !!! |
|---|
| 217 | |
|---|
| 218 | Btus = anuga.Dirichlet_boundary([20.0, 0, 0]) # Outflow water at 10.0 |
|---|
| 219 | Btds = anuga.Dirichlet_boundary([19.0, 0, 0]) # Outflow water at 9.0 |
|---|
| 220 | domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br}) |
|---|
| 221 | |
|---|
| 222 | |
|---|
| 223 | |
|---|
| 224 | #------------------------------------------------------------------------------ |
|---|
| 225 | # Evolve system through time |
|---|
| 226 | #------------------------------------------------------------------------------ |
|---|
| 227 | |
|---|
| 228 | for t in domain.evolve(yieldstep = 1, finaltime = 100): |
|---|
| 229 | print domain.timestepping_statistics() |
|---|
| 230 | |
|---|
| 231 | domain.print_operator_timestepping_statistics() |
|---|
| 232 | |
|---|
