source: anuga_work/development/anuga_1d/channel_domain.py @ 6042

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

Found bug associated with changing max_speed_array to an integer

File size: 17.2 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 domain import *
47Generic_Domain = Domain #Rename
48
49#Shallow water domain
50class Domain(Generic_Domain):
51
52    def __init__(self, coordinates, boundary = None, tagged_elements = None):
53
54        conserved_quantities = ['stage', 'xmomentum']
55        evolved_quantities = ['stage', 'xmomentum', 'elevation', 'height', 'velocity']
56        other_quantities = ['friction']
57        Generic_Domain.__init__(self,
58                                coordinates = coordinates,
59                                boundary = boundary,
60                                conserved_quantities = conserved_quantities,
61                                evolved_quantities = evolved_quantities,
62                                other_quantities = other_quantities,
63                                tagged_elements = tagged_elements)
64       
65        from config import minimum_allowed_height, g, h0
66        self.minimum_allowed_height = minimum_allowed_height
67        self.g = g
68        self.h0 = h0
69
70        #forcing terms not included in 1d domain ?WHy?
71        self.forcing_terms.append(gravity)
72        #self.forcing_terms.append(manning_friction)
73        #print "\nI have Removed forcing terms line 64 1dsw"
74
75       
76        #Stored output
77        self.store = True
78        self.format = 'sww'
79        self.smooth = True
80
81       
82        #Reduction operation for get_vertex_values
83        from util import mean
84        self.reduction = mean
85        #self.reduction = min  #Looks better near steep slopes
86
87        self.set_quantities_to_be_stored(['stage','xmomentum'])
88
89        self.__doc__ = 'shallow_water_domain'
90
91        self.check_integrity()
92
93
94    def check_integrity(self):
95
96        #Check that we are solving the shallow water wave equation
97
98        msg = 'First conserved quantity must be "stage"'
99        assert self.conserved_quantities[0] == 'stage', msg
100        msg = 'Second conserved quantity must be "xmomentum"'
101        assert self.conserved_quantities[1] == 'xmomentum', msg
102
103        msg = 'First evolved quantity must be "stage"'
104        assert self.evolved_quantities[0] == 'stage', msg
105        msg = 'Second evolved quantity must be "xmomentum"'
106        assert self.evolved_quantities[1] == 'xmomentum', msg
107        msg = 'Third evolved quantity must be "elevation"'
108        assert self.evolved_quantities[2] == 'elevation', msg
109        msg = 'Fourth evolved quantity must be "height"'
110        assert self.evolved_quantities[3] == 'height', msg
111        msg = 'Fifth evolved quantity must be "velocity"'
112        assert self.evolved_quantities[4] == 'velocity', msg
113
114        Generic_Domain.check_integrity(self)
115
116    def compute_fluxes(self):
117        #Call correct module function
118        #(either from this module or C-extension)
119        compute_fluxes_vel(self)
120
121    def distribute_to_vertices_and_edges(self):
122        #Call correct module function
123        #(either from this module or C-extension)
124        distribute_to_vertices_and_edges_limit_w_u(self)
125       
126
127#=============== End of Channel Domain ===============================
128#-----------------------------------
129# Compute fluxes interface
130#-----------------------------------
131def compute_fluxes(domain):
132    """
133    Python version of compute fluxes (local_compute_fluxes)
134    is available in test_shallow_water_vel_domain.py
135    """
136 
137
138    from Numeric import zeros, Float
139    import sys
140
141       
142    timestep = float(sys.maxint)
143
144    stage = domain.quantities['stage']
145    xmom = domain.quantities['xmomentum']
146    bed = domain.quantities['elevation']
147
148
149    from comp_flux_vel_ext import compute_fluxes_ext
150
151    domain.flux_timestep = compute_fluxes_ext(timestep,domain,stage,xmom,bed)
152
153
154#-----------------------------------
155# Compute flux definition with vel
156#-----------------------------------
157def compute_fluxes_vel(domain):
158    from Numeric import zeros, Float
159    import sys
160
161       
162    timestep = float(sys.maxint)
163
164    stage    = domain.quantities['stage']
165    xmom     = domain.quantities['xmomentum']
166    bed      = domain.quantities['elevation']
167    height   = domain.quantities['height']
168    velocity = domain.quantities['velocity']
169
170
171    from comp_flux_vel_ext import compute_fluxes_vel_ext
172
173    domain.flux_timestep = compute_fluxes_vel_ext(timestep,domain,stage,xmom,bed,height,velocity)
174
175
176#--------------------------------------------------------------------------
177def distribute_to_vertices_and_edges_limit_w_u(domain):
178    """Distribution from centroids to vertices specific to the
179    shallow water wave
180    equation.
181
182    It will ensure that h (w-z) is always non-negative even in the
183    presence of steep bed-slopes by taking a weighted average between shallow
184    and deep cases.
185
186    In addition, all conserved quantities get distributed as per either a
187    constant (order==1) or a piecewise linear function (order==2).
188
189    FIXME: more explanation about removal of artificial variability etc
190
191    Precondition:
192      All quantities defined at centroids and bed elevation defined at
193      vertices.
194
195    Postcondition
196      Conserved quantities defined at vertices
197
198    """
199
200    #from config import optimised_gradient_limiter
201
202    #Remove very thin layers of water
203    #protect_against_infinitesimal_and_negative_heights(domain) 
204
205    import sys
206    from Numeric import zeros, Float
207    from config import epsilon, h0
208       
209    N = domain.number_of_elements
210
211    #Shortcuts
212    Stage = domain.quantities['stage']
213    Xmom = domain.quantities['xmomentum']
214    Bed = domain.quantities['elevation']
215    Height = domain.quantities['height']
216    Velocity = domain.quantities['velocity']
217
218    #Arrays   
219    w_C   = Stage.centroid_values   
220    uh_C  = Xmom.centroid_values   
221    z_C   = Bed.centroid_values
222    h_C   = Height.centroid_values
223    u_C   = Velocity.centroid_values
224       
225    #print id(h_C)
226    for i in range(N):
227        h_C[i] = w_C[i] - z_C[i]                                               
228        if h_C[i] <= 1.0e-12:
229            #print 'h_C[%d]= %15.5e\n' % (i,h_C[i])
230            h_C[i] = 0.0
231            w_C[i] = z_C[i]
232            #uh_C[i] = 0.0
233           
234 #           u_C[i] = 0.0
235 #       else:
236 #           u_C[i] = uh_C[i]/h_C[i]
237               
238    h0 = 1.0e-12   
239    for i in range(len(h_C)):
240        if h_C[i] < 1.0e-12:
241            u_C[i] = 0.0  #Could have been negative
242            h_C[i] = 0.0
243        else:
244            u_C[i]  = uh_C[i]/(h_C[i] + h0/h_C[i])
245            #u_C[i]  = uh_C[i]/h_C[i]
246       
247    for name in [ 'velocity', 'stage' ]:
248        Q = domain.quantities[name]
249        if domain.order == 1:
250            Q.extrapolate_first_order()
251        elif domain.order == 2:
252            Q.extrapolate_second_order()
253        else:
254            raise 'Unknown order'
255
256    w_V  = domain.quantities['stage'].vertex_values                     
257    z_V  = domain.quantities['elevation'].vertex_values
258    h_V  = domain.quantities['height'].vertex_values
259    u_V  = domain.quantities['velocity'].vertex_values         
260    uh_V = domain.quantities['xmomentum'].vertex_values
261
262    h_V[:]    = w_V - z_V
263    for i in range(len(h_C)):
264        for j in range(2):
265            if h_V[i,j] < 0.0 :
266                #print 'h_V[%d,%d] = %f \n' % (i,j,h_V[i,j])                 
267                dh = h_V[i,j]
268                h_V[i,j] = 0.0
269                w_V[i,j] = z_V[i,j]
270                h_V[i,(j+1)%2] = h_V[i,(j+1)%2] + dh
271                w_V[i,(j+1)%2] = w_V[i,(j+1)%2] + dh
272               
273    uh_V[:] = u_V * h_V
274
275   
276    return
277
278#---------------------------------------------------------------------------
279def distribute_to_vertices_and_edges_limit_w_uh(domain):
280    """Distribution from centroids to vertices specific to the
281    shallow water wave equation.
282
283    In addition, all conserved quantities get distributed as per either a
284    constant (order==1) or a piecewise linear function (order==2).
285
286    Precondition:
287      All quantities defined at centroids and bed elevation defined at
288      vertices.
289
290    Postcondition
291      Conserved quantities defined at vertices
292
293    """
294
295    import sys
296    from Numeric import zeros, Float
297    from config import epsilon, h0
298       
299    N = domain.number_of_elements
300
301    #Shortcuts
302    Stage = domain.quantities['stage']
303    Xmom = domain.quantities['xmomentum']
304    Bed = domain.quantities['elevation']
305    Height = domain.quantities['height']
306    Velocity = domain.quantities['velocity']
307
308    #Arrays   
309    w_C   = Stage.centroid_values   
310    uh_C  = Xmom.centroid_values   
311    z_C   = Bed.centroid_values
312    h_C   = Height.centroid_values
313    u_C   = Velocity.centroid_values
314       
315
316    for i in range(N):
317        h_C[i] = w_C[i] - z_C[i]                                               
318        if h_C[i] <= 1.0e-6:
319            #print 'h_C[%d]= %15.5e\n' % (i,h_C[i])
320            h_C[i] = 0.0
321            w_C[i] = z_C[i]
322            uh_C[i] = 0.0
323           
324    for name in [ 'stage', 'xmomentum']:
325        Q = domain.quantities[name]
326        if domain.order == 1:
327            Q.extrapolate_first_order()
328        elif domain.order == 2:
329            Q.extrapolate_second_order()
330        else:
331            raise 'Unknown order'
332
333    w_V  = domain.quantities['stage'].vertex_values                     
334    z_V  = domain.quantities['elevation'].vertex_values
335    h_V  = domain.quantities['height'].vertex_values
336    u_V  = domain.quantities['velocity'].vertex_values         
337    uh_V = domain.quantities['xmomentum'].vertex_values
338
339    h_V[:]    = w_V - z_V
340
341    for i in range(len(h_C)):
342        for j in range(2):
343            if h_V[i,j] < 0.0 :
344                #print 'h_V[%d,%d] = %f \n' % (i,j,h_V[i,j])                 
345                dh = h_V[i,j]
346                h_V[i,j] = 0.0
347                w_V[i,j] = z_V[i,j]
348                h_V[i,(j+1)%2] = h_V[i,(j+1)%2] + dh
349                w_V[i,(j+1)%2] = w_V[i,(j+1)%2] + dh
350                u_V[i,j] = 0.0
351            if h_V[i,j] < h0:
352                u_V[i,j]
353            else:
354                u_V[i,j] = uh_V[i,j]/(h_V[i,j] + h0/h_V[i,j])
355
356
357#--------------------------------------------------------
358#Boundaries - specific to the shallow_water_vel_domain
359#--------------------------------------------------------
360class Reflective_boundary(Boundary):
361    """Reflective boundary returns same conserved quantities as
362    those present in its neighbour volume but reflected.
363
364    This class is specific to the shallow water equation as it
365    works with the momentum quantities assumed to be the second
366    and third conserved quantities.
367    """
368
369    def __init__(self, domain = None):
370        Boundary.__init__(self)
371
372        if domain is None:
373            msg = 'Domain must be specified for reflective boundary'
374            raise msg
375
376        #Handy shorthands
377        self.normals  = domain.normals
378        self.stage    = domain.quantities['stage'].vertex_values
379        self.xmom     = domain.quantities['xmomentum'].vertex_values
380        self.bed      = domain.quantities['elevation'].vertex_values
381        self.height   = domain.quantities['height'].vertex_values
382        self.velocity = domain.quantities['velocity'].vertex_values
383
384        from Numeric import zeros, Float
385        #self.conserved_quantities = zeros(3, Float)
386        self.evolved_quantities = zeros(5, Float)
387
388    def __repr__(self):
389        return 'Reflective_boundary'
390
391
392    def evaluate(self, vol_id, edge_id):
393        """Reflective boundaries reverses the outward momentum
394        of the volume they serve.
395        """
396
397        q = self.evolved_quantities
398        q[0] = self.stage[vol_id, edge_id]
399        q[1] = -self.xmom[vol_id, edge_id]
400        q[2] = self.bed[vol_id, edge_id]
401        q[3] = self.height[vol_id, edge_id]
402        q[4] = -self.velocity[vol_id, edge_id]
403
404        #print "In Reflective q ",q
405
406
407        return q
408
409class Dirichlet_boundary(Boundary):
410    """Dirichlet boundary returns constant values for the
411    conserved quantities
412    """
413
414
415    def __init__(self, evolved_quantities=None):
416        Boundary.__init__(self)
417
418        if evolved_quantities is None:
419            msg = 'Must specify one value for each evolved quantity'
420            raise msg
421
422        from Numeric import array, Float
423        self.evolved_quantities=array(evolved_quantities).astype(Float)
424
425    def __repr__(self):
426        return 'Dirichlet boundary (%s)' %self.evolved_quantities
427
428    def evaluate(self, vol_id=None, edge_id=None):
429        return self.evolved_quantities
430
431
432#----------------------------
433#Standard forcing terms:
434#---------------------------
435def gravity(domain):
436    """Apply gravitational pull in the presence of bed slope
437    """
438
439    from util import gradient
440    from Numeric import zeros, Float, array, sum
441
442
443
444    Stage     = domain.quantities['stage']
445    Xmom      = domain.quantities['xmomentum']
446    Elevation = domain.quantities['elevation']
447    Height    = domain.quantities['height']
448
449    xmom_ud  = Xmom.explicit_update
450    #stage_ud = Stage.explicit_update
451
452
453    #h = Stage.vertex_values - Elevation.vertex_values
454    h = Height.vertex_values
455    b = Elevation.vertex_values
456    w = Stage.vertex_values
457
458    x = domain.get_vertex_coordinates()
459    g = domain.g
460
461    for k in range(domain.number_of_elements):
462        avg_h = 0.5*(h[k,0] + h[k,1])
463
464        #Compute bed slope
465        x0, x1 = x[k,:]
466        b0, b1 = b[k,:]
467        bx = gradient(x0, x1, b0, b1)
468       
469        #Update momentum (explicit update is reset to source values)
470        xmom_ud[k] += -g*bx*avg_h
471        #stage_ud[k] = 0.0
472 
473 
474def manning_friction(domain):
475    """Apply (Manning) friction to water momentum
476    """
477
478    from math import sqrt
479
480    w = domain.quantities['stage'].centroid_values
481    z = domain.quantities['elevation'].centroid_values
482    h = w-z
483
484    uh = domain.quantities['xmomentum'].centroid_values
485    #vh = domain.quantities['ymomentum'].centroid_values
486    eta = domain.quantities['friction'].centroid_values
487
488    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
489    #ymom_update = domain.quantities['ymomentum'].semi_implicit_update
490
491    N = domain.number_of_elements
492    eps = domain.minimum_allowed_height
493    g = domain.g
494
495    for k in range(N):
496        if eta[k] >= eps:
497            if h[k] >= eps:
498                #S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
499                S = -g * eta[k]**2 * uh[k]
500                S /= h[k]**(7.0/3)
501
502                #Update momentum
503                xmom_update[k] += S*uh[k]
504                #ymom_update[k] += S*vh[k]
505
506def linear_friction(domain):
507    """Apply linear friction to water momentum
508
509    Assumes quantity: 'linear_friction' to be present
510    """
511
512    from math import sqrt
513
514    w = domain.quantities['stage'].centroid_values
515    z = domain.quantities['elevation'].centroid_values
516    h = w-z
517
518    uh = domain.quantities['xmomentum'].centroid_values
519    tau = domain.quantities['linear_friction'].centroid_values
520
521    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
522
523    N = domain.number_of_elements
524    eps = domain.minimum_allowed_height
525
526    for k in range(N):
527        if tau[k] >= eps:
528            if h[k] >= eps:
529                S = -tau[k]/h[k]
530
531                #Update momentum
532                xmom_update[k] += S*uh[k]
533
534
535
536def check_forcefield(f):
537    """Check that f is either
538    1: a callable object f(t,x,y), where x and y are vectors
539       and that it returns an array or a list of same length
540       as x and y
541    2: a scalar
542    """
543
544    from Numeric import ones, Float, array
545
546
547    if callable(f):
548        #N = 3
549        N = 2
550        #x = ones(3, Float)
551        #y = ones(3, Float)
552        x = ones(2, Float)
553        #y = ones(2, Float)
554       
555        try:
556            #q = f(1.0, x=x, y=y)
557            q = f(1.0, x=x)
558        except Exception, e:
559            msg = 'Function %s could not be executed:\n%s' %(f, e)
560            #FIXME: Reconsider this semantics
561            raise msg
562
563        try:
564            q = array(q).astype(Float)
565        except:
566            msg = 'Return value from vector function %s could ' %f
567            msg += 'not be converted into a Numeric array of floats.\n'
568            msg += 'Specified function should return either list or array.'
569            raise msg
570
571        #Is this really what we want?
572        msg = 'Return vector from function %s ' %f
573        msg += 'must have same lenght as input vectors'
574        assert len(q) == N, msg
575
576    else:
577        try:
578            f = float(f)
579        except:
580            msg = 'Force field %s must be either a scalar' %f
581            msg += ' or a vector function'
582            raise msg
583    return f
584
Note: See TracBrowser for help on using the repository browser.