source: anuga_work/development/flow_1d/sww_flow/sww_vel_domain.py @ 7830

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

Moving channel code to numpy

File size: 11.0 KB
Line 
1"""Class Domain -
21D interval domains for finite-volume computations of
3the shallow water wave equation.
4
5This module contains a specialisation of class Domain from module domain.py
6consisting of methods specific to the Shallow Water Wave Equation
7
8
9U_t + E_x = S
10
11where
12
13U = [w, uh]
14E = [uh, u^2h + gh^2/2]
15S represents source terms forcing the system
16(e.g. gravity, friction, wind stress, ...)
17
18and _t, _x, _y denote the derivative with respect to t, x and y respectiely.
19
20The quantities are
21
22symbol    variable name    explanation
23x         x                horizontal distance from origin [m]
24z         elevation        elevation of bed on which flow is modelled [m]
25h         height           water height above z [m]
26w         stage            absolute water level, w = z+h [m]
27u                          speed in the x direction [m/s]
28uh        xmomentum        momentum in the x direction [m^2/s]
29
30eta                        mannings friction coefficient [to appear]
31nu                         wind stress coefficient [to appear]
32
33The conserved quantities are w, uh
34
35For details see e.g.
36Christopher Zoppou and Stephen Roberts,
37Catastrophic Collapse of Water Supply Reservoirs in Urban Areas,
38Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999
39
40
41John Jakeman, Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
42Geoscience Australia, 2006
43"""
44
45
46from flow_1d.generic_1d.generic_domain import *
47from sww_boundary_conditions import *
48from sww_forcing_terms import *
49
50#Shallow water domain
51class Domain(Generic_domain):
52
53    def __init__(self, coordinates, boundary = None, tagged_elements = None):
54
55        conserved_quantities = ['stage', 'xmomentum']
56        evolved_quantities = ['stage', 'xmomentum', 'elevation', 'height', 'velocity']
57        other_quantities = ['friction']
58        Generic_domain.__init__(self,
59                                coordinates = coordinates,
60                                boundary = boundary,
61                                conserved_quantities = conserved_quantities,
62                                evolved_quantities = evolved_quantities,
63                                other_quantities = other_quantities,
64                                tagged_elements = tagged_elements)
65       
66        from flow_1d.config import minimum_allowed_height, g, h0
67        self.minimum_allowed_height = minimum_allowed_height
68        self.g = g
69        self.h0 = h0
70
71        #forcing terms not included in 1d domain ?WHy?
72        self.forcing_terms.append(gravity)
73        #self.forcing_terms.append(manning_friction)
74        #print "\nI have Removed forcing terms line 64 1dsw"
75
76       
77        #Stored output
78        self.store = True
79        self.format = 'sww'
80        self.smooth = True
81
82       
83        #Reduction operation for get_vertex_values
84        from flow_1d.utilities.util import mean
85        self.reduction = mean
86        #self.reduction = min  #Looks better near steep slopes
87
88        self.set_quantities_to_be_stored(['stage','xmomentum'])
89
90        self.__doc__ = 'sww_vel_domain'
91
92        self.check_integrity()
93
94
95
96    def check_integrity(self):
97
98        #Check that we are solving the shallow water wave equation
99
100        msg = 'First conserved quantity must be "stage"'
101        assert self.conserved_quantities[0] == 'stage', msg
102        msg = 'Second conserved quantity must be "xmomentum"'
103        assert self.conserved_quantities[1] == 'xmomentum', msg
104
105        msg = 'First evolved quantity must be "stage"'
106        assert self.evolved_quantities[0] == 'stage', msg
107        msg = 'Second evolved quantity must be "xmomentum"'
108        assert self.evolved_quantities[1] == 'xmomentum', msg
109        msg = 'Third evolved quantity must be "elevation"'
110        assert self.evolved_quantities[2] == 'elevation', msg
111        msg = 'Fourth evolved quantity must be "height"'
112        assert self.evolved_quantities[3] == 'height', msg
113        msg = 'Fifth evolved quantity must be "velocity"'
114        assert self.evolved_quantities[4] == 'velocity', msg
115
116        Generic_domain.check_integrity(self)
117
118    def compute_fluxes(self):
119        #Call correct module function
120        #(either from this module or C-extension)
121        compute_fluxes_vel(self)
122
123    def distribute_to_vertices_and_edges(self):
124        #Call correct module function
125        #(either from this module or C-extension)
126        distribute_to_vertices_and_edges_limit_w_u(self)
127       
128
129#=============== End of Shallow Water Domain ===============================
130#-----------------------------------
131# Compute fluxes interface
132#-----------------------------------
133def compute_fluxes(domain):
134    """
135    Python version of compute fluxes (local_compute_fluxes)
136    is available in test_shallow_water_vel_domain.py
137    """
138 
139
140    from numpy import zeros
141    import sys
142
143       
144    timestep = float(sys.maxint)
145
146    stage = domain.quantities['stage']
147    xmom = domain.quantities['xmomentum']
148    bed = domain.quantities['elevation']
149
150
151    from sww_vel_comp_flux_ext import compute_fluxes_ext
152
153    domain.flux_timestep = compute_fluxes_ext(timestep,domain,stage,xmom,bed)
154
155
156#-----------------------------------
157# Compute flux definition with vel
158#-----------------------------------
159def compute_fluxes_vel(domain):
160    from Numeric import zeros, Float
161    import sys
162
163       
164    timestep = float(sys.maxint)
165
166    stage    = domain.quantities['stage']
167    xmom     = domain.quantities['xmomentum']
168    bed      = domain.quantities['elevation']
169    height   = domain.quantities['height']
170    velocity = domain.quantities['velocity']
171
172
173    from sww_vel_comp_flux_ext import compute_fluxes_vel_ext
174
175    domain.flux_timestep = compute_fluxes_vel_ext(timestep,domain,stage,xmom,bed,height,velocity)
176
177
178
179#--------------------------------------------------------------------------
180def distribute_to_vertices_and_edges_limit_w_u(domain):
181    """Distribution from centroids to vertices specific to the
182    shallow water wave
183    equation.
184
185    It will ensure that h (w-z) is always non-negative even in the
186    presence of steep bed-slopes by taking a weighted average between shallow
187    and deep cases.
188
189    In addition, all conserved quantities get distributed as per either a
190    constant (order==1) or a piecewise linear function (order==2).
191
192    FIXME: more explanation about removal of artificial variability etc
193
194    Precondition:
195      All quantities defined at centroids and bed elevation defined at
196      vertices.
197
198    Postcondition
199      Conserved quantities defined at vertices
200
201    """
202
203    #from config import optimised_gradient_limiter
204
205    #Remove very thin layers of water
206    #protect_against_infinitesimal_and_negative_heights(domain) 
207
208    import sys
209    from numpy import zeros
210    from config import epsilon, h0
211       
212    N = domain.number_of_elements
213
214    #Shortcuts
215    Stage = domain.quantities['stage']
216    Xmom = domain.quantities['xmomentum']
217    Bed = domain.quantities['elevation']
218    Height = domain.quantities['height']
219    Velocity = domain.quantities['velocity']
220
221    #Arrays   
222    w_C   = Stage.centroid_values   
223    uh_C  = Xmom.centroid_values   
224    z_C   = Bed.centroid_values
225    h_C   = Height.centroid_values
226    u_C   = Velocity.centroid_values
227       
228    #print id(h_C)
229##     for i in range(N):
230##      h_C[i] = w_C[i] - z_C[i]                                               
231##         if h_C[i] <= 1.0e-12:
232##          #print 'h_C[%d]= %15.5e\n' % (i,h_C[i])
233##          h_C[i] = 0.0
234##          w_C[i] = z_C[i]
235##          #uh_C[i] = 0.0
236           
237##  #           u_C[i] = 0.0
238##  #       else:
239##  #           u_C[i] = uh_C[i]/h_C[i]
240               
241    h0 = 1.0e-12   
242    for i in range(N):
243        h_C[i] = w_C[i] - z_C[i]                                               
244        if h_C[i] < 1.0e-12:
245            u_C[i] = 0.0  #Could have been negative
246            h_C[i] = 0.0
247            w_C[i] = z_C[i]
248        else:
249            #u_C[i]  = uh_C[i]/(h_C[i] + h0/h_C[i])
250            u_C[i]  = uh_C[i]/h_C[i]
251       
252    for name in [ 'velocity', 'stage' ]:
253        Q = domain.quantities[name]
254        if domain.order == 1:
255            Q.extrapolate_first_order()
256        elif domain.order == 2:
257            Q.extrapolate_second_order()
258        else:
259            raise 'Unknown order'
260
261    w_V  = domain.quantities['stage'].vertex_values                     
262    z_V  = domain.quantities['elevation'].vertex_values
263    h_V  = domain.quantities['height'].vertex_values
264    u_V  = domain.quantities['velocity'].vertex_values         
265    uh_V = domain.quantities['xmomentum'].vertex_values
266
267    #print w_V
268    #print z_V
269   
270    h_V[:,:]    = w_V - z_V
271    for i in range(len(h_C)):
272        for j in range(2):
273            if h_V[i,j] < 0.0 :
274                #print 'h_V[%d,%d] = %f \n' % (i,j,h_V[i,j])                 
275                dh = h_V[i,j]
276                h_V[i,j] = 0.0
277                w_V[i,j] = z_V[i,j]
278                h_V[i,(j+1)%2] = h_V[i,(j+1)%2] + dh
279                w_V[i,(j+1)%2] = w_V[i,(j+1)%2] + dh
280               
281    uh_V[:,:] = u_V * h_V
282
283   
284    return
285
286#---------------------------------------------------------------------------
287def distribute_to_vertices_and_edges_limit_w_uh(domain):
288    """Distribution from centroids to vertices specific to the
289    shallow water wave equation.
290
291    In addition, all conserved quantities get distributed as per either a
292    constant (order==1) or a piecewise linear function (order==2).
293
294    Precondition:
295      All quantities defined at centroids and bed elevation defined at
296      vertices.
297
298    Postcondition
299      Conserved quantities defined at vertices
300
301    """
302
303    import sys
304    from Numeric import zeros, Float
305    from config import epsilon, h0
306       
307    N = domain.number_of_elements
308
309    #Shortcuts
310    Stage = domain.quantities['stage']
311    Xmom = domain.quantities['xmomentum']
312    Bed = domain.quantities['elevation']
313    Height = domain.quantities['height']
314    Velocity = domain.quantities['velocity']
315
316    #Arrays   
317    w_C   = Stage.centroid_values   
318    uh_C  = Xmom.centroid_values   
319    z_C   = Bed.centroid_values
320    h_C   = Height.centroid_values
321    u_C   = Velocity.centroid_values
322       
323
324    for i in range(N):
325        h_C[i] = w_C[i] - z_C[i]                                               
326        if h_C[i] <= 1.0e-6:
327            #print 'h_C[%d]= %15.5e\n' % (i,h_C[i])
328            h_C[i] = 0.0
329            w_C[i] = z_C[i]
330            uh_C[i] = 0.0
331           
332    for name in [ 'stage', 'xmomentum']:
333        Q = domain.quantities[name]
334        if domain.order == 1:
335            Q.extrapolate_first_order()
336        elif domain.order == 2:
337            Q.extrapolate_second_order()
338        else:
339            raise 'Unknown order'
340
341    w_V  = domain.quantities['stage'].vertex_values                     
342    z_V  = domain.quantities['elevation'].vertex_values
343    h_V  = domain.quantities['height'].vertex_values
344    u_V  = domain.quantities['velocity'].vertex_values         
345    uh_V = domain.quantities['xmomentum'].vertex_values
346
347    h_V[:,:]    = w_V - z_V
348
349    for i in range(len(h_C)):
350        for j in range(2):
351            if h_V[i,j] < 0.0 :
352                #print 'h_V[%d,%d] = %f \n' % (i,j,h_V[i,j])                 
353                dh = h_V[i,j]
354                h_V[i,j] = 0.0
355                w_V[i,j] = z_V[i,j]
356                h_V[i,(j+1)%2] = h_V[i,(j+1)%2] + dh
357                w_V[i,(j+1)%2] = w_V[i,(j+1)%2] + dh
358                u_V[i,j] = 0.0
359            if h_V[i,j] < h0:
360                u_V[i,j]
361            else:
362                u_V[i,j] = uh_V[i,j]/(h_V[i,j] + h0/h_V[i,j])
363
364
Note: See TracBrowser for help on using the repository browser.