source: trunk/anuga_core/source/anuga/structures/test_Outlet_Ctrl.py @ 8020

Last change on this file since 8020 was 8020, checked in by steve, 14 years ago

Fixed bug with log_to_file on boyd_box_operator

  • Property svn:executable set to *
File size: 9.8 KB
RevLine 
[7997]1""" Testing CULVERT (Changing from Horizontal Abstraction to Vertical Abstraction
2
3This example includes a Model Topography that shows a TYPICAL Headwall Configuration
4
5The aim is to change the Culvert Routine to Model more precisely the abstraction
6from a vertical face.
7
8The inflow must include the impact of Approach velocity.
9Similarly the Outflow has MOMENTUM Not just Up welling as in the Horizontal Style
10abstraction
11
12"""
13print 'Starting.... Importing Modules...'
14#------------------------------------------------------------------------------
15# Import necessary modules
16#------------------------------------------------------------------------------
17import anuga
18
19from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
20from anuga.shallow_water.shallow_water_domain import Domain
21from 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
[7998]27from anuga.structures.boyd_pipe_operator import Boyd_pipe_operator
[8020]28from anuga.structures.boyd_box_operator import Boyd_box_operator
[7997]29#from anuga.culvert_flows.culvert_routines import weir_orifice_channel_culvert_model
30from math import pi,pow,sqrt
31
32import numpy as num
33
34
35
36
37#------------------------------------------------------------------------------
38# Setup computational domain
39#------------------------------------------------------------------------------
40print 'Setting up domain'
41
42length = 120. #x-Dir
43width = 200.  #y-dir
44
45dx = 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
50points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
51                                                    len1=length, len2=width)
52domain = Domain(points, vertices, boundary)   
53domain.set_name('Test_Outlet_Ctrl')                 # Output name
54domain.set_default_order(2)
55domain.H0 = 0.01
56domain.tight_slope_limiters = 1
57
58print 'Size', len(domain)
59
60#------------------------------------------------------------------------------
61# Setup initial conditions
62#------------------------------------------------------------------------------
63
64def 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
127print 'Setting Quantities....'
128domain.set_quantity('elevation', topography)  # Use function for elevation
129domain.set_quantity('friction', 0.01)         # Constant friction
130domain.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#------------------------------------------------------------------------------
[8020]143print 'Defining Structures'
[7997]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#                            height=10.0,
153#                            apron=5.0,
154#                            verbose=False)
155
156
157#------------------------------------------------------------------------------
158# Setup culverts
159#------------------------------------------------------------------------------
160
161culverts = []
[8018]162number_of_culverts = 2
[7997]163for 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 = [40.0, y]
167    ep1 = [50.0, y]
[7998]168    losses = {'inlet':0.5, 'outlet':1, 'bend':0, 'grate':0, 'pier': 0, 'other': 0}
[8020]169#    culverts.append(Boyd_pipe_operator(domain,
170#                            end_point0=ep0,
171#                            end_point1=ep1,
172#                            losses=losses,
173#                            diameter=1.5, #culvert_width, #3.658,
174#                            apron=6.0,
175#                            use_momentum_jet=True,
176#                            use_velocity_head=True,
177#                            manning=0.013,
178#                            logging=True,
179#                            label='culvert',
180#                            verbose=False))
181
182
183#    Boyd_pipe_operator(domain,
184#                            end_point0=ep0,
185#                            end_point1=ep1,
186#                            losses=losses,
187#                            diameter=1.5, #culvert_width, #3.658,
188#                            apron=6.0,
189#                            use_momentum_jet=True,
190#                            use_velocity_head=True,
191#                            manning=0.013,
192#                            logging=True,
193#                            label='pipe_culvert',
194#                            verbose=False)
195
196    Boyd_box_operator(domain,
[7997]197                            end_point0=ep0,
198                            end_point1=ep1,
[7998]199                            losses=losses,
[8020]200                            width=culvert_width,
201                            height=10.0,
[7997]202                            apron=6.0,
203                            use_momentum_jet=True,
204                            use_velocity_head=True,
205                            manning=0.013,
[8018]206                            logging=True,
[8020]207                            label='box_culvert',
208                            verbose=False)
[7997]209
[7998]210                       
[7997]211
[7998]212#losses = {'inlet':1, 'outlet':1, 'bend':1, 'grate':1, 'pier': 1, 'other': 1}
213#culvert2 = Culvert_operator(domain,
214                            #end_point0=[40.0, 62.5],
215                            #end_point1=[50.0, 62.5],
216                            #losses,
217                            #width=25.0,
218                            #height=10.0,
219                            #apron=5.0,
220                            #manning=0.013,
221                            #verbose=False)
[7997]222
223
224
225#------------------------------------------------------------------------------
226# Setup boundary conditions
227#------------------------------------------------------------------------------
228print 'Setting Boundary Conditions'
229Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
230Bi = anuga.Dirichlet_boundary([0.0, 0.0, 0.0])          # Inflow based on Flow Depth and Approaching Momentum !!!
231
232Btus = anuga.Dirichlet_boundary([20.0, 0, 0])           # Outflow water at 10.0
233Btds = anuga.Dirichlet_boundary([19.0, 0, 0])           # Outflow water at 9.0
234domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br})
235
236
[8018]237
[7997]238#------------------------------------------------------------------------------
239# Evolve system through time
240#------------------------------------------------------------------------------
241
242for t in domain.evolve(yieldstep = 1, finaltime = 100):
243    print domain.timestepping_statistics()
244
[8018]245    domain.print_operator_timestepping_statistics()
246
Note: See TracBrowser for help on using the repository browser.