source: anuga_work/development/flow_1d/channel_flow/channel_domain_Ab.py @ 7832

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

Putting 1d stuff in one folder

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