source: branches/numpy/anuga/culvert_flows/Test_Culvert_Flat_Water_Lev.py @ 6533

Last change on this file since 6533 was 6533, checked in by rwilson, 15 years ago

Revert back to 6481, prior to auto-merge of trunk and numpy branch.

File size: 5.9 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#------------------------------------------------------------------------------
17from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
18
19from anuga.shallow_water import Domain, Reflective_boundary,\
20     Dirichlet_boundary,\
21     Transmissive_boundary, Time_boundary
22
23from anuga.culvert_flows.culvert_class import Culvert_flow
24from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
25     
26from math import pi,pow,sqrt
27
28import numpy as num
29
30
31#------------------------------------------------------------------------------
32# Setup computational domain
33#------------------------------------------------------------------------------
34print 'Setting up domain'
35
36length = 40.
37width = 5.
38
39dx = dy = 1           # Resolution: Length of subdivisions on both axes
40#dx = dy = .5           # Resolution: Length of subdivisions on both axes
41#dx = dy = .5           # Resolution: Length of subdivisions on both axes
42#dx = dy = .1           # Resolution: Length of subdivisions on both axes
43
44points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
45                                               len1=length, len2=width)
46domain = Domain(points, vertices, boundary)   
47domain.set_name('Test_Culv_Flat_WL')                 # Output name
48domain.set_default_order(2)
49domain.H0 = 0.01
50domain.tight_slope_limiters = 1
51
52print 'Size', len(domain)
53
54#------------------------------------------------------------------------------
55# Setup initial conditions
56#------------------------------------------------------------------------------
57
58def topography(x, y):
59    """Set up a weir
60   
61    A culvert will connect either side
62    """
63    # General Slope of Topography
64    z=-x/1000
65   
66    #       NOW Add bits and Pieces to topography
67    N = len(x)
68    for i in range(N):
69
70       # Sloping Embankment Across Channel
71        if 5.0 < x[i] < 10.1:
72            if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: # Cut Out Segment for Culvert FACE
73               z[i]=z[i]
74            else:
75               z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
76        if 10.0 < x[i] < 12.1:
77           z[i] +=  2.5              # Flat Crest of Embankment
78        if 12.0 < x[i] < 14.5:
79            if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5: # Cut Out Segment for Culvert FACE
80               z[i]=z[i]
81            else:
82               z[i] +=  2.5-1.0*(x[i] -12.0)       # Sloping D/S Face
83                           
84       
85               
86    return z
87
88print 'Setting Quantities....'
89domain.set_quantity('elevation', topography)  # Use function for elevation
90domain.set_quantity('friction', 0.01)         # Constant friction
91domain.set_quantity('stage',
92                    expression='elevation')   # Dry initial condition
93
94
95
96
97#------------------------------------------------------------------------------
98# Setup specialised forcing terms
99#------------------------------------------------------------------------------
100
101#------------------------------------------------------------------------------
102# Setup CULVERT INLETS and OUTLETS in Current Topography
103#------------------------------------------------------------------------------
104print 'DEFINING any Structures if Required'
105
106#  DEFINE CULVERT INLET AND OUTLETS
107
108
109culvert_rating = Culvert_flow(domain,
110                       culvert_description_filename='example_rating_curve.csv',
111                       end_point0=[9.0, 2.5], 
112                       end_point1=[13.0, 2.5],
113                       verbose=True)
114
115
116culvert_energy = Culvert_flow(domain,
117                       label='Culvert No. 1',
118                       description='This culvert is a test unit 1.2m Wide by 0.75m High',   
119                       end_point0=[9.0, 2.5], 
120                       end_point1=[13.0, 2.5],
121                       width=1.20,height=0.75,
122                       culvert_routine=boyd_generalised_culvert_model,       
123                       number_of_barrels=1,
124                       update_interval=2,
125                       verbose=True)
126
127domain.forcing_terms.append(culvert_energy)
128
129#------------------------------------------------------------------------------
130# Setup boundary conditions
131#------------------------------------------------------------------------------
132print 'Setting Boundary Conditions'
133Bi = Dirichlet_boundary([0.0, 0.0, 0.0])          # Inflow based on Flow Depth and Approaching Momentum !!!
134Br = Reflective_boundary(domain)              # Solid reflective wall
135Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
136Btus = Time_boundary(domain, lambda t: [0.0+ 1.25*(1+num.sin(2*pi*(t-4)/10)), 0.0, 0.0])
137Btds = Time_boundary(domain, lambda t: [0.0+ 0.75*(1+num.sin(2*pi*(t-4)/20)), 0.0, 0.0])
138domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br})
139
140
141#------------------------------------------------------------------------------
142# Evolve system through time
143#------------------------------------------------------------------------------
144
145#for t in domain.evolve(yieldstep = 1, finaltime = 25):
146#    print domain.timestepping_statistics()
147   
148
149   
150   
151#import sys; sys.exit()
152# Profiling code
153import time
154t0 = time.time()
155   
156s = 'for t in domain.evolve(yieldstep = 1, finaltime = 25): domain.write_time()'
157
158import profile, pstats
159FN = 'profile.dat'
160
161profile.run(s, FN)
162   
163print 'That took %.2f seconds' %(time.time()-t0)
164
165S = pstats.Stats(FN)
166#S.sort_stats('time').print_stats(20)
167s = S.sort_stats('cumulative').print_stats(30)
168
169print s
Note: See TracBrowser for help on using the repository browser.