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

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

Added new tests

  • Property svn:executable set to *
File size: 8.0 KB
Line 
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
26from anuga.culvert_flows.culvert_class import Culvert_flow
27from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
28#from anuga.culvert_flows.culvert_routines import weir_orifice_channel_culvert_model
29from math import pi,pow,sqrt
30
31import numpy as num
32
33
34#------------------------------------------------------------------------------
35# Setup computational domain
36#------------------------------------------------------------------------------
37print 'Setting up domain'
38
39length = 200. #x-Dir
40width = 200.  #y-dir
41
42dx = dy = 2.5          # Resolution: Length of subdivisions on both axes
43#dx = dy = .5           # Resolution: Length of subdivisions on both axes
44#dx = dy = .5           # Resolution: Length of subdivisions on both axes
45#dx = dy = .1           # Resolution: Length of subdivisions on both axes
46
47points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
48                                               len1=length, len2=width)
49domain = Domain(points, vertices, boundary)   
50domain.set_name('Test_Open_slot_WIDE_BRIDGE')                 # Output name
51domain.set_default_order(2)
52domain.H0 = 0.01
53domain.tight_slope_limiters = 1
54
55print 'Size', len(domain)
56
57#------------------------------------------------------------------------------
58# Setup initial conditions
59#------------------------------------------------------------------------------
60
61def topography(x, y):
62    """Set up a weir
63   
64    A culvert will connect either side
65    """
66    # General Slope of Topography
67    z=10.0-x/100.0  # % Longitudinal Slope
68   
69    #       NOW Add bits and Pieces to topography
70    bank_hgt=10.0
71    bridge_width = 50.0
72    bank_width = 10.0
73   
74    us_apron_skew = 1.0 # 1.0 = 1 Length: 1 Width, 2.0 = 2 Length : 1 Width
75    us_start_x = 10.0
76    top_start_y = 50.0
77    us_slope = 3.0  #Horiz : Vertic
78    ds_slope = 3.0
79    ds_apron_skew = 1.0 # 1.0 = 1 Length: 1 Width, 2.0 = 2 Length : 1 Width
80    centre_line_y= top_start_y+bridge_width/2.0
81
82    # CALCULATE PARAMETERS TO FORM THE EMBANKMENT
83    us_slope_length = bank_hgt*us_slope
84    us_end_x =us_start_x + us_slope_length
85    us_toe_start_y =top_start_y - us_slope_length / us_apron_skew
86    us_toe_end_y = top_start_y + bridge_width + us_slope_length / us_apron_skew
87
88    top_end_y = top_start_y + bridge_width
89    ds_slope_length = bank_hgt*ds_slope
90    ds_start_x = us_end_x + bank_width
91    ds_end_x = ds_start_x + ds_slope_length
92
93    ds_toe_start_y =top_start_y - ds_slope_length / ds_apron_skew
94    ds_toe_end_y = top_start_y + bridge_width + ds_slope_length / ds_apron_skew
95
96
97    N = len(x)
98    for i in range(N):
99
100       # Sloping Embankment Across Channel
101        if us_start_x < x[i] < us_end_x +0.1:   # For UPSLOPE on the Upstream FACE
102        #if 5.0 < x[i] < 10.1: # For a Range of X, and over a Range of Y based on X adjust Z
103            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:
104                        #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
105               z[i]=z[i] # Flat Apron
106            else:
107               z[i] += z[i] + (x[i] - us_start_x)/us_slope    # Set elevation for Sloping Segment  U/S Face
108        if us_end_x < x[i] < ds_start_x + 0.1: # FOR The TOP of BANK Segment
109           if top_start_y < y[i] < top_end_y:
110               z[i]=z[i] # Flat Apron
111           else:
112               z[i] +=  z[i]+bank_hgt        # Flat Crest of Embankment
113        if ds_start_x < x[i] < ds_end_x: # DOWN SDLOPE Segment on Downstream face
114            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
115               z[i]=z[i] # Flat Apron
116            else:
117               z[i] += z[i]+bank_hgt-(x[i] -ds_start_x)/ds_slope       # Sloping D/S Face
118           
119       
120
121    return z
122
123print 'Setting Quantities....'
124domain.set_quantity('elevation', topography)  # Use function for elevation
125domain.set_quantity('friction', 0.01)         # Constant friction
126domain.set_quantity('stage',
127                    expression='elevation')   # Dry initial condition
128
129
130
131
132#------------------------------------------------------------------------------
133# Setup specialised forcing terms
134#------------------------------------------------------------------------------
135
136#------------------------------------------------------------------------------
137# Setup CULVERT INLETS and OUTLETS in Current Topography
138#------------------------------------------------------------------------------
139print 'DEFINING any Structures if Required'
140
141#  DEFINE CULVERT INLET AND OUTLETS
142
143"""
144culvert_rating = Culvert_flow(domain,
145                       culvert_description_filename='example_rating_curve.csv',
146                       end_point0=[40.0, 75.0],
147                       end_point1=[50.0, 75.0],
148                       verbose=True)
149
150
151culvert_energy = Culvert_flow(domain,
152                       label='Culvert No. 1',
153                       description='This culvert is a test unit 1.2m Wide by 0.75m High',   
154                       end_point0=[40.0, 75.0],
155                       end_point1=[50.0, 75.0],
156                       width=50.0,height=5.0,
157                       culvert_routine=boyd_generalised_culvert_model,       
158                       number_of_barrels=1,
159                       #update_interval=0.25,
160                       log_file=True,
161                       discharge_hydrograph=True,
162                       use_velocity_head=False,
163                       verbose=True)
164
165domain.forcing_terms.append(culvert_energy)
166"""
167
168#------------------------------------------------------------------------------
169# Setup boundary conditions
170#------------------------------------------------------------------------------
171print 'Setting Boundary Conditions'
172Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
173Bi = anuga.Dirichlet_boundary([0.0, 0.0, 0.0])          # Inflow based on Flow Depth and Approaching Momentum !!!
174
175Bo = anuga.Dirichlet_boundary([-5.0, 0, 0])           # Outflow water at -5.0
176Bd = anuga.Dirichlet_boundary([0,0,0])                # Outflow water at 0.0
177
178#Btus = Time_boundary(domain, lambda t: [0.0+ 1.025*(1+num.sin(2*pi*(t-4)/10)), 0.0, 0.0])
179#Btds = Time_boundary(domain, lambda t: [0.0+ 0.0075*(1+num.sin(2*pi*(t-4)/20)), 0.0, 0.0])
180
181Btus = anuga.Dirichlet_boundary([18.0, 0, 0])           # Outflow water at 5.0
182Btds = anuga.Dirichlet_boundary([0.0, 0, 0])           # Outflow water at 1.0
183domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br})
184
185
186#------------------------------------------------------------------------------
187# Evolve system through time
188#------------------------------------------------------------------------------
189
190for t in domain.evolve(yieldstep = 1, finaltime = 100):
191    print domain.timestepping_statistics()
Note: See TracBrowser for help on using the repository browser.