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

Last change on this file since 8018 was 8018, checked in by steve, 13 years ago

Added logging to file for fractional step operators. There needs to an member function log_timestepping_statistics() implemented for each operator.

  • Property svn:executable set to *
File size: 8.6 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
26#from anuga.culvert_flows.culvert_class import Culvert_flow
27from anuga.structures.boyd_pipe_operator import Boyd_pipe_operator
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
36#------------------------------------------------------------------------------
37# Setup computational domain
38#------------------------------------------------------------------------------
39print 'Setting up domain'
40
41length = 120. #x-Dir
42width = 200.  #y-dir
43
44dx = dy = 2.0          # Resolution: Length of subdivisions on both axes
45#dx = dy = .5           # Resolution: Length of subdivisions on both axes
46#dx = dy = .5           # Resolution: Length of subdivisions on both axes
47#dx = dy = .1           # Resolution: Length of subdivisions on both axes
48
49points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
50                                                    len1=length, len2=width)
51domain = Domain(points, vertices, boundary)   
52domain.set_name('Test_Outlet_Ctrl')                 # Output name
53domain.set_default_order(2)
54domain.H0 = 0.01
55domain.tight_slope_limiters = 1
56
57print 'Size', len(domain)
58
59#------------------------------------------------------------------------------
60# Setup initial conditions
61#------------------------------------------------------------------------------
62
63def topography(x, y):
64    """Set up a weir
65   
66    A culvert will connect either side
67    """
68    # General Slope of Topography
69    z=10.0-x/100.0  # % Longitudinal Slope
70   
71    #       NOW Add bits and Pieces to topography
72    bank_hgt=10.0
73    bridge_width = 50.0
74    bank_width = 10.0
75   
76    us_apron_skew = 1.0 # 1.0 = 1 Length: 1 Width, 2.0 = 2 Length : 1 Width
77    us_start_x = 10.0
78    top_start_y = 50.0
79    us_slope = 3.0  #Horiz : Vertic
80    ds_slope = 3.0
81    ds_apron_skew = 1.0 # 1.0 = 1 Length: 1 Width, 2.0 = 2 Length : 1 Width
82    centre_line_y= top_start_y+bridge_width/2.0
83
84    # CALCULATE PARAMETERS TO FORM THE EMBANKMENT
85    us_slope_length = bank_hgt*us_slope
86    us_end_x =us_start_x + us_slope_length
87    us_toe_start_y =top_start_y - us_slope_length / us_apron_skew
88    us_toe_end_y = top_start_y + bridge_width + us_slope_length / us_apron_skew
89
90    top_end_y = top_start_y + bridge_width
91    ds_slope_length = bank_hgt*ds_slope
92    ds_start_x = us_end_x + bank_width
93    ds_end_x = ds_start_x + ds_slope_length
94
95    ds_toe_start_y =top_start_y - ds_slope_length / ds_apron_skew
96    ds_toe_end_y = top_start_y + bridge_width + ds_slope_length / ds_apron_skew
97
98
99    N = len(x)
100    for i in range(N):
101
102       # Sloping Embankment Across Channel
103        if us_start_x < x[i] < us_end_x +0.1:   # For UPSLOPE on the Upstream FACE
104        #if 5.0 < x[i] < 10.1: # For a Range of X, and over a Range of Y based on X adjust Z
105            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:
106                #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
107                 z[i]=z[i] # Flat Apron
108                #z[i] += z[i] + (x[i] - us_start_x)/us_slope
109                #pass
110            else:
111               z[i] += z[i] + (x[i] - us_start_x)/us_slope    # Sloping Segment  U/S Face
112        if us_end_x < x[i] < ds_start_x + 0.1:
113           z[i] +=  z[i]+bank_hgt        # Flat Crest of Embankment
114        if ds_start_x < x[i] < ds_end_x: # DOWN SDLOPE Segment on Downstream face
115            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
116                 z[i]=z[i] # Flat Apron
117                #z[i] += z[i]+bank_hgt-(x[i] -ds_start_x)/ds_slope
118                #pass
119            else:
120               z[i] += z[i]+bank_hgt-(x[i] -ds_start_x)/ds_slope       # Sloping D/S Face
121           
122       
123
124    return z
125
126print 'Setting Quantities....'
127domain.set_quantity('elevation', topography)  # Use function for elevation
128domain.set_quantity('friction', 0.01)         # Constant friction
129domain.set_quantity('stage',
130                    expression='elevation')   # Dry initial condition
131
132
133
134
135#------------------------------------------------------------------------------
136# Setup specialised forcing terms
137#------------------------------------------------------------------------------
138
139#------------------------------------------------------------------------------
140# Setup CULVERT INLETS and OUTLETS in Current Topography
141#------------------------------------------------------------------------------
142print 'DEFINING any Structures if Required'
143
144#  DEFINE CULVERT INLET AND OUTLETS
145
146
147#culvert0 = Culvert_operator(domain,
148#                            end_point0=[40.0, 75.0],
149#                            end_point1=[50.0, 75.0],
150#                            width=50.0,
151#                            height=10.0,
152#                            apron=5.0,
153#                            verbose=False)
154
155
156#------------------------------------------------------------------------------
157# Setup culverts
158#------------------------------------------------------------------------------
159
160culverts = []
161number_of_culverts = 2
162for i in range(number_of_culverts):
163    culvert_width = 50.0/number_of_culverts
164    y = 100-i*culvert_width - culvert_width/2.0
165    ep0 = [40.0, y]
166    ep1 = [50.0, y]
167    losses = {'inlet':0.5, 'outlet':1, 'bend':0, 'grate':0, 'pier': 0, 'other': 0}
168    culverts.append(Boyd_pipe_operator(domain,
169                            end_point0=ep0,
170                            end_point1=ep1,
171                            losses=losses,
172                            diameter=1.5, #culvert_width, #3.658,
173                            apron=6.0,
174                            use_momentum_jet=True,
175                            use_velocity_head=True,
176                            manning=0.013,
177                            logging=True,
178                            verbose=False))
179
180                       
181
182#losses = {'inlet':1, 'outlet':1, 'bend':1, 'grate':1, 'pier': 1, 'other': 1}
183#culvert2 = Culvert_operator(domain,
184                            #end_point0=[40.0, 62.5],
185                            #end_point1=[50.0, 62.5],
186                            #losses,
187                            #width=25.0,
188                            #height=10.0,
189                            #apron=5.0,
190                            #manning=0.013,
191                            #verbose=False)
192
193
194
195#------------------------------------------------------------------------------
196# Setup boundary conditions
197#------------------------------------------------------------------------------
198print 'Setting Boundary Conditions'
199Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
200Bi = anuga.Dirichlet_boundary([0.0, 0.0, 0.0])          # Inflow based on Flow Depth and Approaching Momentum !!!
201
202Btus = anuga.Dirichlet_boundary([20.0, 0, 0])           # Outflow water at 10.0
203Btds = anuga.Dirichlet_boundary([19.0, 0, 0])           # Outflow water at 9.0
204domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br})
205
206
207
208#------------------------------------------------------------------------------
209# Evolve system through time
210#------------------------------------------------------------------------------
211
212for t in domain.evolve(yieldstep = 1, finaltime = 100):
213    print domain.timestepping_statistics()
214
215    domain.print_operator_timestepping_statistics()
216
Note: See TracBrowser for help on using the repository browser.