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

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

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

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