source: anuga_work/development/anuga_1d/backup.py @ 6453

Last change on this file since 6453 was 6453, checked in by steve, 15 years ago

Added Padarn's (2008/2009 summer student) files

File size: 16.1 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
8This particular modification of the Domain class implements the ability to
9vary the width of the 1D channel that the water flows in. As a result the
10conserved variables are different than previous implementations and so are the
11equations.
12
13U_t + E_x = S
14
15where
16------------!!!! NOTE THIS NEEDS UPDATING !!!!------------------
17U = [A, Q]
18E = [Q, Q^2/A + gh^2/2]
19S represents source terms forcing the system
20(e.g. gravity, friction, wind stress, ...)
21
22and _t, _x, _y denote the derivative with respect to t, x and y respectiely.
23
24The quantities are
25
26symbol    variable name    explanation
27x         x                horizontal distance from origin [m]
28z         elevation        elevation of bed on which flow is modelled [m]
29h         height           water height above z [m]
30w         stage            absolute water level, w = z+h [m]
31u                          speed in the x direction [m/s]
32uh        xmomentum        momentum in the x direction [m^2/s]
33
34eta                        mannings friction coefficient [to appear]
35nu                         wind stress coefficient [to appear]
36
37The conserved quantities are w, uh
38--------------------------------------------------------------------------
39For details see e.g.
40Christopher Zoppou and Stephen Roberts,
41Catastrophic Collapse of Water Supply Reservoirs in Urban Areas,
42Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999
43
44
45John Jakeman, Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou,
46Padarn Wilson, Geoscience Australia, 2008
47"""
48
49
50from domain import *
51Generic_Domain = Domain #Rename
52
53#Shallow water domain
54class Domain(Generic_Domain):
55
56    def __init__(self, coordinates, boundary = None, tagged_elements = None):
57
58        conserved_quantities = ['area', 'discharge']
59        evolved_quantities = ['area', 'discharge', 'elevation', 'height', 'velocity','width','stage']
60        other_quantities = ['friction']
61        Generic_Domain.__init__(self,
62                                coordinates = coordinates,
63                                boundary = boundary,
64                                conserved_quantities = conserved_quantities,
65                                evolved_quantities = evolved_quantities,
66                                other_quantities = other_quantities,
67                                tagged_elements = tagged_elements)
68       
69        from config import minimum_allowed_height, g, h0
70        self.minimum_allowed_height = minimum_allowed_height
71        self.g = g
72        self.h0 = h0
73
74     
75        self.forcing_terms.append(gravity)
76        self.forcing_terms.append(boundary_stress)
77        #self.forcing_terms.append(manning_friction)
78       
79
80       
81        #Stored output
82        self.store = True
83        self.format = 'sww'
84        self.smooth = True
85
86       
87        #Reduction operation for get_vertex_values
88        from util import mean
89        self.reduction = mean
90        #self.reduction = min  #Looks better near steep slopes
91
92        self.set_quantities_to_be_stored(['area','discharge'])
93
94        self.__doc__ = 'channel_domain_Ab'
95
96        self.check_integrity()
97
98
99    def check_integrity(self):
100
101        #Check that we are solving the shallow water wave equation
102
103        msg = 'First conserved quantity must be "area"'
104        assert self.conserved_quantities[0] == 'area', msg
105        msg = 'Second conserved quantity must be "discharge"'
106        assert self.conserved_quantities[1] == 'discharge', msg
107
108        msg = 'First evolved quantity must be "area"'
109        assert self.evolved_quantities[0] == 'area', msg
110        msg = 'Second evolved quantity must be "discharge"'
111        assert self.evolved_quantities[1] == 'discharge', msg
112        msg = 'Third evolved quantity must be "elevation"'
113        assert self.evolved_quantities[2] == 'elevation', msg
114        msg = 'Fourth evolved quantity must be "height"'
115        assert self.evolved_quantities[3] == 'height', msg
116        msg = 'Fifth evolved quantity must be "velocity"'
117        assert self.evolved_quantities[4] == 'velocity', msg
118        msg = 'Fifth evolved quantity must be "width"'
119        assert self.evolved_quantities[5] == 'width', msg
120
121        Generic_Domain.check_integrity(self)
122
123    def compute_fluxes(self):
124        #Call correct module function
125        #(either from this module or C-extension)
126        compute_fluxes_channel(self)
127
128    def distribute_to_vertices_and_edges(self):
129        #Call correct module function
130        #(either from this module or C-extension)
131        distribute_to_vertices_and_edges_limit_a_d(self)
132       
133
134#=============== End of Channel Domain ===============================
135
136#-----------------------------------
137# Compute flux definition with channel
138#-----------------------------------
139def compute_fluxes_channel(domain):
140    from Numeric import zeros, Float
141    import sys
142
143       
144    timestep = float(sys.maxint)
145
146    area    = domain.quantities['area']
147    discharge     = domain.quantities['discharge']
148    bed      = domain.quantities['elevation']
149    height   = domain.quantities['height']
150    velocity = domain.quantities['velocity']
151    width    = domain.quantities['width']
152
153   
154 
155    from channel_domain_ext import compute_fluxes_channel_ext
156    domain.flux_timestep = compute_fluxes_channel_ext(timestep,domain,area,discharge,bed,height,velocity,width)
157
158#-----------------------------------------------------------------------
159# Distribute to verticies with stage reconstructed and then extrapolated
160#-----------------------------------------------------------------------
161def distribute_to_vertices_and_edges_limit_a_d(domain):
162   
163    #Remove very thin layers of water
164    #protect_against_infinitesimal_and_negative_heights(domain) 
165
166    import sys
167    from Numeric import zeros, Float
168    from config import epsilon, h0
169       
170    N = domain.number_of_elements
171
172    #Shortcuts
173    Area = domain.quantities['area']
174    Discharge = domain.quantities['discharge']
175    Bed = domain.quantities['elevation']
176    Height = domain.quantities['height']
177    Velocity = domain.quantities['velocity']
178    Width = domain.quantities['width']
179    Stage = domain.quantities['stage']
180
181    #Arrays   
182    a_C   = Area.centroid_values   
183    d_C  = Discharge.centroid_values   
184    z_C   = Bed.centroid_values
185    h_C   = Height.centroid_values
186    u_C   = Velocity.centroid_values
187    b_C   = Width.centroid_values
188    w_C   = Stage.centroid_values
189
190    h0 = 1.0e-12       
191    #print id(h_C)
192    for i in range(N):
193                                                       
194        if a_C[i] <= h0:
195            a_C[i] = 0.0
196            h_C[i] = 0.0
197            d_C[i] = 0.0
198            u_C[i] = 0.0
199            w_C[i] = z_C[i]
200           
201           
202
203        else:
204
205            if b_C[i]<=h0:
206                a_C[i] = 0.0
207                h_C[i] = 0.0
208                d_C[i] = 0.0
209                u_C[i] = 0.0
210                w_C[i] = z_C[i]
211
212            else:
213                h_C[i] = a_C[i]/(b_C[i]+h0/b_C[i])
214                w_C[i] = h_C[i]+z_C[i]
215                u_C[i] = d_C[i]/(a_C[i]+h0/a_C[i])
216               
217   
218       
219    for name in ['discharge','elevation','width','stage','height','area','velocity']:
220        Q = domain.quantities[name]
221        if domain.order == 1:
222            Q.extrapolate_first_order()
223        elif domain.order == 2:
224            Q.extrapolate_second_order()
225        else:
226            raise 'Unknown order'
227    a_V  = domain.quantities['area'].vertex_values
228    w_V  = domain.quantities['stage'].vertex_values                     
229    z_V  = domain.quantities['elevation'].vertex_values
230    h_V  = domain.quantities['height'].vertex_values
231    u_V  = domain.quantities['velocity'].vertex_values         
232    d_V = domain.quantities['discharge'].vertex_values
233    b_V = domain.quantities['width'].vertex_values
234
235   
236##     for i in range(len(h_C)):
237##         for j in range(2):
238##             if b_V[i,j] < h0 :
239##                 a_V[i,j]=0
240##                 h_V[i,j]=0
241##                 d_V[i,j]=0
242##                 u_V[i,j]=0
243##             else:                   
244##                 h_V[i,j] = w_V[i,j]-z_V[i,j]
245##                 if h_V[i,j]<h0:
246##                     h_V[i,j]=0
247##                     w_V[i,j]=z_V[i,j]
248##                 a_V[i,j] = b_V[i,j]*h_V[i,j]
249##                 if d_V[i,j]<h0 or a_V[i,j]<h0:
250##                     d_V[i,j]=0
251##                     u_V[i,j]=0
252##                 else:
253##                     u_V[i,j]=d_V[i,j]/(a_V[i,j]+h0/a_V[i,j])
254           
255               
256               
257   
258
259   
260    return
261
262
263#--------------------------------------------------------
264#Boundaries - specific to the shallow_water_vel_domain
265#--------------------------------------------------------
266class Reflective_boundary(Boundary):
267    """Reflective boundary returns same conserved quantities as
268    those present in its neighbour volume but reflected.
269
270    This class is specific to the shallow water equation as it
271    works with the momentum quantities assumed to be the second
272    and third conserved quantities.
273    """
274
275    def __init__(self, domain = None):
276        Boundary.__init__(self)
277
278        if domain is None:
279            msg = 'Domain must be specified for reflective boundary'
280            raise msg
281
282        #Handy shorthands
283        self.normals  = domain.normals
284        self.area    = domain.quantities['area'].vertex_values
285        self.discharge     = domain.quantities['discharge'].vertex_values
286        self.bed      = domain.quantities['elevation'].vertex_values
287        self.height   = domain.quantities['height'].vertex_values
288        self.velocity = domain.quantities['velocity'].vertex_values
289        self.width    = domain.quantities['width'].vertex_values
290        self.stage    = domain.quantities['stage'].vertex_values
291
292        from Numeric import zeros, Float
293        #self.conserved_quantities = zeros(3, Float)
294        self.evolved_quantities = zeros(7, Float)
295
296    def __repr__(self):
297        return 'Reflective_boundary'
298
299
300    def evaluate(self, vol_id, edge_id):
301        """Reflective boundaries reverses the outward momentum
302        of the volume they serve.
303        """
304
305        q = self.evolved_quantities
306        q[0] = self.area[vol_id, edge_id]
307        q[1] = -self.discharge[vol_id, edge_id]
308        q[2] = self.bed[vol_id, edge_id]
309        q[3] = self.height[vol_id, edge_id]
310        q[4] = -self.velocity[vol_id, edge_id]
311        q[5] = self.width[vol_id,edge_id]
312        q[6] = self.stage[vol_id,edge_id]
313
314        #print "In Reflective q ",q
315
316
317        return q
318
319class Dirichlet_boundary(Boundary):
320    """Dirichlet boundary returns constant values for the
321    conserved quantities
322if k>5 and k<15:
323            print discharge_ud[k],-g*zx*avg_h*avg_b
324        discharge_ud[k] +=-g*zx*avg_h*avg_b    """
325
326
327    def __init__(self, evolved_quantities=None):
328        Boundary.__init__(self)
329
330        if evolved_quantities is None:
331            msg = 'Must specify one value for each evolved quantity'
332            raise msg
333
334        from Numeric import array, Float
335        self.evolved_quantities=array(evolved_quantities).astype(Float)
336
337    def __repr__(self):
338        return 'Dirichlet boundary (%s)' %self.evolved_quantities
339
340    def evaluate(self, vol_id=None, edge_id=None):
341        return self.evolved_quantities
342
343
344#----------------------------
345#Standard forcing terms:
346#---------------------------
347def gravity(domain):
348    """Apply gravitational pull in the presence of bed slope
349    """
350
351    from util import gradient
352    from Numeric import zeros, Float, array, sum
353
354
355
356    Area     = domain.quantities['area']
357    Discharge  = domain.quantities['discharge']
358    Elevation = domain.quantities['elevation']
359    Height    = domain.quantities['height']
360    Width     = domain.quantities['width']
361
362    discharge_ud  = Discharge.explicit_update
363
364
365       
366    h = Height.vertex_values
367    b = Width.vertex_values
368    a = Area.vertex_values
369    z = Elevation.vertex_values
370
371    x = domain.get_vertex_coordinates()
372    g = domain.g
373 
374    for k in range(domain.number_of_elements):
375        avg_h = 0.5*(h[k,0] + h[k,1])
376        avg_b = 0.5*(b[k,0] + b[k,1])
377       
378        #Compute bed slope
379        x0, x1 = x[k,:]
380        z0, z1 = z[k,:]
381        zx = gradient(x0, x1, z0, z1)
382        #Update momentum (explicit update is reset to source values)
383        discharge_ud[k]+= -g*zx*avg_h*avg_b
384       
385     
386
387def boundary_stress(domain):
388
389
390    from util import gradient
391    from Numeric import zeros, Float, array, sum
392
393
394
395    Area     = domain.quantities['area']
396    Discharge  = domain.quantities['discharge']
397    Elevation = domain.quantities['elevation']
398    Height    = domain.quantities['height']
399    Width     = domain.quantities['width']
400
401    discharge_ud  = Discharge.explicit_update
402
403
404       
405    h = Height.vertex_values
406    b = Width.vertex_values
407    a = Area.vertex_values
408    z = Elevation.vertex_values
409
410    x = domain.get_vertex_coordinates()
411    g = domain.g
412
413    for k in range(domain.number_of_elements):
414        avg_h = 0.5*(h[k,0] + h[k,1])
415       
416
417        #Compute bed slope
418        x0, x1 = x[k,:]
419        b0, b1 = b[k,:]
420        bx = gradient(x0, x1, b0, b1)
421       
422        #Update momentum (explicit update is reset to source values)
423        discharge_ud[k] += 0.5*g*bx*avg_h*avg_h
424        #stage_ud[k] = 0.0
425 
426 
427def manning_friction(domain):
428    """Apply (Manning) friction to water momentum
429    """
430
431    from math import sqrt
432
433    w = domain.quantities['stage'].centroid_values
434    z = domain.quantities['elevation'].centroid_values
435    h = w-z
436
437    uh = domain.quantities['xmomentum'].centroid_values
438    #vh = domain.quantities['ymomentum'].centroid_values
439    eta = domain.quantities['friction'].centroid_values
440
441    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
442    #ymom_update = domain.quantities['ymomentum'].semi_implicit_update
443
444    N = domain.number_of_elements
445    eps = domain.minimum_allowed_height
446    g = domain.g
447
448    for k in range(N):
449        if eta[k] >= eps:
450            if h[k] >= eps:
451                #S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
452                S = -g * eta[k]**2 * uh[k]
453                S /= h[k]**(7.0/3)
454
455                #Update momentum
456                xmom_update[k] += S*uh[k]
457                #ymom_update[k] += S*vh[k]
458
459def linear_friction(domain):
460    """Apply linear friction to water momentum
461
462    Assumes quantity: 'linear_friction' to be present
463    """
464
465    from math import sqrt
466
467    w = domain.quantities['stage'].centroid_values
468    z = domain.quantities['elevation'].centroid_values
469    h = w-z
470
471    uh = domain.quantities['xmomentum'].centroid_values
472    tau = domain.quantities['linear_friction'].centroid_values
473
474    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
475
476    N = domain.number_of_elements
477    eps = domain.minimum_allowed_height
478
479    for k in range(N):
480        if tau[k] >= eps:
481            if h[k] >= eps:
482                S = -tau[k]/h[k]
483
484                #Update momentum
485                xmom_update[k] += S*uh[k]
486
487
488
489def check_forcefield(f):
490    """Check that f is either
491    1: a callable object f(t,x,y), where x and y are vectors
492       and that it returns an array or a list of same length
493       as x and y
494    2: a scalar
495    """
496
497    from Numeric import ones, Float, array
498
499
500    if callable(f):
501        #N = 3
502        N = 2
503        #x = ones(3, Float)
504        #y = ones(3, Float)
505        x = ones(2, Float)
506        #y = ones(2, Float)
507       
508        try:
509            #q = f(1.0, x=x, y=y)
510            q = f(1.0, x=x)
511        except Exception, e:
512            msg = 'Function %s could not be executed:\n%s' %(f, e)
513            #FIXME: Reconsider this semantics
514            raise msg
515
516        try:
517            q = array(q).astype(Float)
518        except:
519            msg = 'Return value from vector function %s could ' %f
520            msg += 'not be converted into a Numeric array of floats.\n'
521            msg += 'Specified function should return either list or array.'
522            raise msg
523
524        #Is this really what we want?
525        msg = 'Return vector from function %s ' %f
526        msg += 'must have same lenght as input vectors'
527        assert len(q) == N, msg
528
529    else:
530        try:
531            f = float(f)
532        except:
533            msg = 'Force field %s must be either a scalar' %f
534            msg += ' or a vector function'
535            raise msg
536    return f
537
Note: See TracBrowser for help on using the repository browser.