# source:anuga_work/development/Rudy_vandrie_2007/culvert_flow.py@4697

Last change on this file since 4697 was 4697, checked in by ole, 16 years ago

Culvert flow prototype for Rudy

File size: 8.9 KB
Line
1"""Example of shallow water wave equation.
2
3Specific methods pertaining to the 2D shallow water equation
4are imported from shallow_water
5for use with the generic finite volume framework
6
7Conserved quantities are h, uh and vh stored as elements 0, 1 and 2 in the
8numerical vector named conserved_quantities.
9"""
10
11#------------------------------------------------------------------------------
12# Import necessary modules
13#------------------------------------------------------------------------------
14from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
15from anuga.shallow_water import Domain
16from anuga.shallow_water.shallow_water_domain import Reflective_boundary
17from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
18
19
20#------------------------------------------------------------------------------
21# Setup computational domain
22#------------------------------------------------------------------------------
23length = 40.
24width = 5.
25
26#dx = dy = 1           # Resolution: Length of subdivisions on both axes
27#dx = dy = .5           # Resolution: Length of subdivisions on both axes
28dx = dy = .1           # Resolution: Length of subdivisions on both axes
29
30points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
31                                               len1=length, len2=width)
32domain = Domain(points, vertices, boundary)
33domain.set_name('culvert_example')                 # Output name
34print 'Size', len(domain)
35
36#------------------------------------------------------------------------------
37# Setup initial conditions
38#------------------------------------------------------------------------------
39
40def topography(x, y):
41    """Set up a weir
42
43    A culvert will connect either side
44    """
45
46    z = -x/10
47
48    N = len(x)
49    for i in range(N):
50
51        # Weir from 10 to 12 m
52        if 10 < x[i] < 12:
53            z[i] += 0.4 - 0.05*y[i]
54
55        # Constriction
56        #if 27 < x[i] < 29 and y[i] > 3:
57        #    z[i] += 2
58
59        # Pole
60        #if (x[i] - 34)**2 + (y[i] - 2)**2 < 0.4**2:
61        #    z[i] += 2
62
63    return z
64
65
66domain.set_quantity('elevation', topography)  # Use function for elevation
67domain.set_quantity('friction', 0.01)         # Constant friction
68domain.set_quantity('stage',
69                    expression='elevation')   # Dry initial condition
70
71
72
73
74#------------------------------------------------------------------------------
75# Setup specialised forcing terms
76#------------------------------------------------------------------------------
77
78class Inflow:
79    """Class Inflow - general 'rain and drain' forcing term.
80
81    Useful for implementing flows in and out of the domain.
82
84
85    center [m]: Coordinates at center of flow point
86    radius [m]: Size of circular area
87    flow [m/s]: Rate of change of quantity over the specified area.
88                This parameter can be either a constant or a function of time.
89                Positive values indicate inflow,
90                negative values indicate outflow.
91
92    Examples
93
94    Inflow((0.7, 0.4), 0.07, -0.2) # Constant drain at 0.2 m/s.
95                                   # This corresponds to a flow of
96                                   # 0.07**2*pi*0.2 = 0.00314 m^3/s
97
98    Inflow((0.5, 0.5), 0.001, lambda t: min(4*t, 5)) # Tap turning up to
99                                                     # a maximum inflow of
100                                                     # 5 m/s over the
101                                                     # specified area
102    """
103
104    # FIXME (OLE): Add a polygon as an alternative.
105    # FIXME (OLE): Generalise to all quantities
106    # FIXME (OLE):
107
108    def __init__(self,
109                 domain,
111                 flow=0.0,
112                 quantity_name = 'stage'):
113
114        from anuga.utilities.numerical_tools import ensure_numeric
115
116        if center is not None and radius is not None:
117            assert len(center) == 2
118        else:
119            msg = 'Both center and radius must be specified'
120            raise Exception, msg
121
122        self.domain = domain
123        self.center = ensure_numeric(center)
125        self.flow = flow
126
127        self.quantity = domain.quantities[quantity_name].explicit_update
128
129        # Determine indices in flow area
130        N = len(domain)
131        self.indices = []
132        coordinates = domain.get_centroid_coordinates()
133        for k in range(N):
134            x, y = coordinates[k,:] # Centroid
136                self.indices.append(k)
137
138
139    def __call__(self, domain):
140
141        # Update inflow
142        if callable(self.flow):
143            flow = self.flow(domain.get_time())
144        else:
145            flow = self.flow
146
147        for k in self.indices:
148            self.quantity[k] += flow
149
150
151class Culvert_flow:
152    """Culvert flow - transfer water from one hole to another
153    """
154
155    def __init__(self,
156                 domain,
159                 transfer_flow=None):
160
161        from Numeric import sqrt, sum
162
163        self.openings = []
164        self.openings.append(Inflow(domain,
165                                    center=center0,
167                                    flow=None))
168
169        self.openings.append(Inflow(domain,
170                                    center=center1,
172                                    flow=None))
173
174        # Compute distance between opening centers. I assume there are only two.
175        self.distance = sqrt(sum( (self.openings[0].center - self.openings[1].center)**2 ))
176        print 'Distance between openings is', self.distance
177
178
179    def __call__(self, domain):
180        from anuga.utilities.numerical_tools import mean
181
182
183        # Get average water depths at each opening
184        for opening in self.openings:
185            stage = domain.quantities['stage'].get_values(location='centroids',
186                                                          indices=opening.indices)
187            elevation = domain.quantities['elevation'].get_values(location='centroids',
188                                                              indices=opening.indices)
189
190            # Store current avg depth with opening object
191            opening.depth = mean(stage-elevation)
192            # print 'Depth', opening.depth
193
194
195        # Here's a simple transfer function hardwired into this Culvert class:
196        # Let flow rate depend on difference in depth
197        # Flow is instantaneous for now but one might use
198        # the precomputed distance between the two openings
199        # to somehow introduce a delay.
200
201        delta_h = self.openings[0].depth - self.openings[1].depth
202        flow_rate = 0.5*9.81*delta_h**2                          # Just an example of flow rate
204
205        if delta_h > 0:
206            # Opening 0 has the greatest depth.
207            # Flow will go from opening 0 to opening 1
208
209            self.openings[0].flow = -flow_rate
210            self.openings[1].flow = flow_rate
211        else:
212            # Opening 1 has the greatest depth.
213            # Flow will go from opening 1 to opening 0
214
215            self.openings[0].flow = flow_rate
216            self.openings[1].flow = -flow_rate
217
218
219        # Execute flow term for each opening
220        for opening in self.openings:
221            opening(domain)
222
223
224
225#sink = Inflow(domain, center=[0.7, 0.4], radius=0.0707, flow = 0.0)
226#tap = Inflow(domain, center=(0.5, 0.5), radius=0.0316,
227#            flow=lambda t: min(4*t, 5)) # Tap turning up
228#domain.forcing_terms.append(tap)
229#domain.forcing_terms.append(sink)
230
231culvert = Culvert_flow(domain,
234                       transfer_flow=None) # Hardwired for now
235
236domain.forcing_terms.append(culvert)
237
238#------------------------------------------------------------------------------
239# Setup boundary conditions
240#------------------------------------------------------------------------------
241Bi = Dirichlet_boundary([0.4, 0, 0])          # Inflow
242Br = Reflective_boundary(domain)              # Solid reflective wall
243Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
244
245domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
246
247
248#------------------------------------------------------------------------------
249# Evolve system through time
250#------------------------------------------------------------------------------
251for t in domain.evolve(yieldstep = 0.01, finaltime = 15):
252    domain.write_time()
253
254    #if domain.get_time() >= 4 and tap.flow != 0.0:
255    #    print 'Turning tap off'
256    #    tap.flow = 0.0
257
258    #if domain.get_time() >= 3 and sink.flow < 0.0:
259    #    print 'Turning drain on'
260    #    sink.flow = -0.8
261
Note: See TracBrowser for help on using the repository browser.