source: anuga_work/development/2010-projects/anuga_1d/sww/sww_vel_domain.py @ 7860

Last change on this file since 7860 was 7860, checked in by steve, 12 years ago

Continuing to numpy the for loops

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 anuga_1d.base.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 anuga_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 anuga_1d.base.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 anuga_1d.sww_flow.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 anuga_1d.sww.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
204    #Remove very thin layers of water
205    #protect_against_infinitesimal_and_negative_heights(domain) 
206
207    import sys
208    from numpy import zeros
209    from anuga_1d.config import epsilon, h0
210       
211    N = domain.number_of_elements
212
213    #Shortcuts
214    Stage = domain.quantities['stage']
215    Xmom = domain.quantities['xmomentum']
216    Bed = domain.quantities['elevation']
217    Height = domain.quantities['height']
218    Velocity = domain.quantities['velocity']
219
220    #Arrays   
221    w_C   = Stage.centroid_values   
222    uh_C  = Xmom.centroid_values   
223    z_C   = Bed.centroid_values
224    h_C   = Height.centroid_values
225    u_C   = Velocity.centroid_values
226       
227    #print id(h_C)
228##     for i in range(N):
229##      h_C[i] = w_C[i] - z_C[i]                                               
230##         if h_C[i] <= 1.0e-12:
231##          #print 'h_C[%d]= %15.5e\n' % (i,h_C[i])
232##          h_C[i] = 0.0
233##          w_C[i] = z_C[i]
234##          #uh_C[i] = 0.0
235           
236##  #           u_C[i] = 0.0
237##  #       else:
238##  #           u_C[i] = uh_C[i]/h_C[i]
239               
240    h0 = 1.0e-12   
241    for i in range(N):
242        h_C[i] = w_C[i] - z_C[i]                                               
243        if h_C[i] < 1.0e-12:
244            u_C[i] = 0.0  #Could have been negative
245            h_C[i] = 0.0
246            w_C[i] = z_C[i]
247        else:
248            u_C[i]  = uh_C[i]/(h_C[i] + h0/h_C[i])
249            #u_C[i]  = uh_C[i]/h_C[i]
250       
251    for name in [ 'velocity', 'stage' ]:
252        Q = domain.quantities[name]
253        if domain.order == 1:
254            Q.extrapolate_first_order()
255        elif domain.order == 2:
256            Q.extrapolate_second_order()
257        else:
258            raise 'Unknown order'
259
260    w_V  = domain.quantities['stage'].vertex_values                     
261    z_V  = domain.quantities['elevation'].vertex_values
262    h_V  = domain.quantities['height'].vertex_values
263    u_V  = domain.quantities['velocity'].vertex_values         
264    uh_V = domain.quantities['xmomentum'].vertex_values
265
266    #print w_V
267    #print z_V
268   
269    h_V[:,:]    = w_V - z_V
270    for i in range(len(h_C)):
271        for j in range(2):
272            if h_V[i,j] < 0.0 :
273                #print 'h_V[%d,%d] = %f \n' % (i,j,h_V[i,j])                 
274                dh = h_V[i,j]
275                h_V[i,j] = 0.0
276                w_V[i,j] = z_V[i,j]
277                h_V[i,(j+1)%2] = h_V[i,(j+1)%2] + dh
278                w_V[i,(j+1)%2] = w_V[i,(j+1)%2] + dh
279               
280    uh_V[:,:] = u_V * h_V
281
282   
283    return
284
285#---------------------------------------------------------------------------
286def distribute_to_vertices_and_edges_limit_w_uh(domain):
287    """Distribution from centroids to vertices specific to the
288    shallow water wave equation.
289
290    In addition, all conserved quantities get distributed as per either a
291    constant (order==1) or a piecewise linear function (order==2).
292
293    Precondition:
294      All quantities defined at centroids and bed elevation defined at
295      vertices.
296
297    Postcondition
298      Conserved quantities defined at vertices
299
300    """
301
302    import sys
303    from numpy import zeros
304    from anuga_1d.config import epsilon, h0
305       
306    N = domain.number_of_elements
307
308    #Shortcuts
309    Stage = domain.quantities['stage']
310    Xmom = domain.quantities['xmomentum']
311    Bed = domain.quantities['elevation']
312    Height = domain.quantities['height']
313    Velocity = domain.quantities['velocity']
314
315    #Arrays   
316    w_C   = Stage.centroid_values   
317    uh_C  = Xmom.centroid_values   
318    z_C   = Bed.centroid_values
319    h_C   = Height.centroid_values
320    u_C   = Velocity.centroid_values
321       
322
323    for i in range(N):
324        h_C[i] = w_C[i] - z_C[i]                                               
325        if h_C[i] <= 1.0e-6:
326            #print 'h_C[%d]= %15.5e\n' % (i,h_C[i])
327            h_C[i] = 0.0
328            w_C[i] = z_C[i]
329            uh_C[i] = 0.0
330           
331    for name in [ 'stage', 'xmomentum']:
332        Q = domain.quantities[name]
333        if domain.order == 1:
334            Q.extrapolate_first_order()
335        elif domain.order == 2:
336            Q.extrapolate_second_order()
337        else:
338            raise 'Unknown order'
339
340    w_V  = domain.quantities['stage'].vertex_values                     
341    z_V  = domain.quantities['elevation'].vertex_values
342    h_V  = domain.quantities['height'].vertex_values
343    u_V  = domain.quantities['velocity'].vertex_values         
344    uh_V = domain.quantities['xmomentum'].vertex_values
345
346    h_V[:,:]    = w_V - z_V
347
348    for i in range(len(h_C)):
349        for j in range(2):
350            if h_V[i,j] < 0.0 :
351                #print 'h_V[%d,%d] = %f \n' % (i,j,h_V[i,j])                 
352                dh = h_V[i,j]
353                h_V[i,j] = 0.0
354                w_V[i,j] = z_V[i,j]
355                h_V[i,(j+1)%2] = h_V[i,(j+1)%2] + dh
356                w_V[i,(j+1)%2] = w_V[i,(j+1)%2] + dh
357                u_V[i,j] = 0.0
358            if h_V[i,j] < h0:
359                u_V[i,j]
360            else:
361                u_V[i,j] = uh_V[i,j]/(h_V[i,j] + h0/h_V[i,j])
362
363
Note: See TracBrowser for help on using the repository browser.