source: anuga_work/development/flow_1d/channel_flow/channel_domain2.py @ 7827

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

Putting 1d stuff in one folder

File size: 26.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 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 = ['area', 'discharge']
55        evolved_quantities = ['area', 'discharge', 'elevation', 'height', 'stage']
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(['area','discharge'])
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_channel(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_a_d(self)
125       
126
127#=============== End of Shallow Water 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    domain.flux_timestep = .1
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#  Compute fluxes channel
177#----------------------------------
178def compute_fluxes_channel(domain):
179    from Numeric import zeros, Float
180    import sys
181   
182    timestep = float(sys.maxint)
183
184    area      = domain.quantities['area']
185    discharge = domain.quantities['discharge']
186    bed       = domain.quantities['elevation']
187    height     = domain.quantities['height']
188    stage     = domain.quantities['stage']
189
190    #from channel_domain_ext import compute_fluxes_channel_ext
191
192   # domain.flux_timestep = compute_fluxes_channel_ext(timestep,domain,area,discharge,bed,height,stage)
193   # domain.quantities['area'].explicit_update=ones(410,Float)
194    domain.flux_timestep = .1
195    #print area.vertex_values[0],area.vertex_values[1],area.vertex_values[2]
196    for i in range(len(domain.coordinates)-1):
197        fluxesl= channel_flux_func(domain,i-1)
198        fluxesr= channel_flux_func(domain,i)
199        #print domain.areas[i]
200        #print -1*(fluxesl[1]-fluxesr[1])
201        #print height.centroid_values[200],height.vertex_values[200]
202        #print fluxesl[0]
203        print -1*(fluxesl[0]-fluxesr[0])*domain.areas[i]
204        area.explicit_update[i]=-1*(fluxesl[0]-fluxesr[0])*domain.areas[i]
205        discharge.explicit_update[i]=-(fluxesl[1]-fluxesr[1])*domain.areas[i]
206       
207def channel_flux_func(domain, i):
208
209    area      = domain.quantities['area']
210    discharge = domain.quantities['discharge']
211    bed       = domain.quantities['elevation']
212    height     = domain.quantities['height']
213    stage     = domain.quantities['stage']
214
215    flux_left0=0
216    flux_left1=0
217    flux_right0=0
218    flux_right1=0
219    flux0=0
220    flux1=0
221  # Crude numerical flux calculation
222    from math import sqrt
223
224    if i==0 or i==(len(domain.coordinates)-1):
225        flux_left0=0
226        flux_left1=0
227        flux_right0=0
228        flux_right1=0
229    else:
230        g=9.8
231        a_left = area.centroid_values[i-1]
232        d_left = discharge.centroid_values[i-1]
233        z_left = bed.centroid_values[i-1]
234        h_left = height.centroid_values[i-1]
235        w_left = stage.centroid_values[i-1]
236
237        a_right = area.centroid_values[i]
238        d_right = discharge.centroid_values[i]
239        z_right = bed.centroid_values[i]
240        h_right = height.centroid_values[i]
241        w_right = stage.centroid_values[i]
242
243        z=(z_left+z_right)/2
244        #hbarr=0.5*(height.vertex_values[i][1]+height.vertex_values[i][0])
245        #hbarl=0.5*(height.vertex_values[i-1][1]+height.vertex_values[i-1][0])
246        hbarr=0
247        hbarl=0
248        #print hbarl,h_left
249        if a_left>1.0e-12:
250            u_left=d_left/a_left
251        else:
252            u_left=0
253        if a_right>1.0e-12:
254            u_right=d_right/a_right
255        else:
256            u_right=0
257
258        ## soundspeed_left = sqrt(g*h_left);
259##         soundspeed_right = sqrt(g*h_right);
260
261##         s_max = max(u_left+soundspeed_left,u_right+soundspeed_right)
262##         if s_max<0.0:
263##             s_max=0
264##         s_min = min(u_left-soundspeed_left,u_right-soundspeed_right)
265##         if s_min>0.0:
266##             s_min=0
267           
268        flux_left0 = d_left
269        flux_right1= d_right
270        if a_left<1.0e-12:
271            flux_left11=0
272        else:
273            flux_left11=d_left*d_left/a_left
274        flux_left1=flux_left11+g*(calculateI(h_left,0)-calculateI(hbarl,0))
275        if a_right<1.0e-12:
276            flux_right11=0
277        else:
278            flux_right11=d_right*d_right/a_right
279        flux_right1=flux_right11+g*(calculateI(h_right,0)-calculateI(hbarr,0))
280        #print g*(calculateI(h_right,0)-calculateI(hbarr,0)), g*(calculateI(h_left,0)-calculateI(hbarl,0))
281
282    flux0 = 0.5*(flux_left0+flux_right0)
283    flux1 = 0.5*(flux_left1+flux_right1)
284    return (flux0,flux1)
285           
286def calculateI(H,z0):
287    return H*H*0.5-z0*H+z0*z0*0.5
288
289   
290#--------------------------------------------------------------------------
291def distribute_to_vertices_and_edges_limit_w_u(domain):
292    """Distribution from centroids to vertices specific to the
293    shallow water wave
294    equation.
295
296    It will ensure that h (w-z) is always non-negative even in the
297    presence of steep bed-slopes by taking a weighted average between shallow
298    and deep cases.
299
300    In addition, all conserved quantities get distributed as per either a
301    constant (order==1) or a piecewise linear function (order==2).
302
303    FIXME: more explanation about removal of artificial variability etc
304
305    Precondition:
306      All quantities defined at centroids and bed elevation defined at
307      vertices.
308
309    Postcondition
310      Conserved quantities defined at vertices
311
312    """
313
314    #from config import optimised_gradient_limiter
315
316    #Remove very thin layers of water
317    #protect_against_infinitesimal_and_negative_heights(domain) 
318
319    import sys
320    from Numeric import zeros, Float
321    from config import epsilon, h0
322       
323    N = domain.number_of_elements
324
325    #Shortcuts
326    Stage = domain.quantities['stage']
327    Xmom = domain.quantities['xmomentum']
328    Bed = domain.quantities['elevation']
329    Height = domain.quantities['height']
330    Velocity = domain.quantities['velocity']
331
332    #Arrays   
333    w_C   = Stage.centroid_values   
334    uh_C  = Xmom.centroid_values   
335    z_C   = Bed.centroid_values
336    h_C   = Height.centroid_values
337    u_C   = Velocity.centroid_values
338       
339    #print id(h_C)
340    for i in range(N):
341        h_C[i] = w_C[i] - z_C[i]                                               
342        if h_C[i] <= 1.0e-12:
343            #print 'h_C[%d]= %15.5e\n' % (i,h_C[i])
344            h_C[i] = 0.0
345            w_C[i] = z_C[i]
346            #uh_C[i] = 0.0
347           
348 #           u_C[i] = 0.0
349 #       else:
350 #           u_C[i] = uh_C[i]/h_C[i]
351               
352    h0 = 1.0e-12   
353    for i in range(len(h_C)):
354        if h_C[i] < 1.0e-12:
355            u_C[i] = 0.0  #Could have been negative
356            h_C[i] = 0.0
357        else:
358            u_C[i]  = uh_C[i]/(h_C[i] + h0/h_C[i])
359            #u_C[i]  = uh_C[i]/h_C[i]
360       
361    for name in [ 'velocity', 'stage' ]:
362        Q = domain.quantities[name]
363        if domain.order == 1:
364            Q.extrapolate_first_order()
365        elif domain.order == 2:
366            Q.extrapolate_second_order()
367        else:
368            raise 'Unknown order'
369
370    w_V  = domain.quantities['stage'].vertex_values                     
371    z_V  = domain.quantities['elevation'].vertex_values
372    h_V  = domain.quantities['height'].vertex_values
373    u_V  = domain.quantities['velocity'].vertex_values         
374    uh_V = domain.quantities['xmomentum'].vertex_values
375
376    h_V[:]    = w_V - z_V
377    for i in range(len(h_C)):
378        for j in range(2):
379            if h_V[i,j] < 0.0 :
380                #print 'h_V[%d,%d] = %f \n' % (i,j,h_V[i,j])                 
381                dh = h_V[i,j]
382                h_V[i,j] = 0.0
383                w_V[i,j] = z_V[i,j]
384                h_V[i,(j+1)%2] = h_V[i,(j+1)%2] + dh
385                w_V[i,(j+1)%2] = w_V[i,(j+1)%2] + dh
386               
387    uh_V[:] = u_V * h_V
388
389   
390    return
391
392#---------------------------------------------------------------------------
393def distribute_to_vertices_and_edges_limit_w_uh(domain):
394    """Distribution from centroids to vertices specific to the
395    shallow water wave equation.
396
397    In addition, all conserved quantities get distributed as per either a
398    constant (order==1) or a piecewise linear function (order==2).
399
400    Precondition:
401      All quantities defined at centroids and bed elevation defined at
402      vertices.
403
404    Postcondition
405      Conserved quantities defined at vertices
406
407    """
408
409    import sys
410    from Numeric import zeros, Float
411    from config import epsilon, h0
412       
413    N = domain.number_of_elements
414
415    #Shortcuts
416    Stage = domain.quantities['stage']
417    Xmom = domain.quantities['xmomentum']
418    Bed = domain.quantities['elevation']
419    Height = domain.quantities['height']
420    Velocity = domain.quantities['velocity']
421
422    #Arrays   
423    w_C   = Stage.centroid_values   
424    uh_C  = Xmom.centroid_values   
425    z_C   = Bed.centroid_values
426    h_C   = Height.centroid_values
427    u_C   = Velocity.centroid_values
428       
429
430    for i in range(N):
431        h_C[i] = w_C[i] - z_C[i]                                               
432        if h_C[i] <= 1.0e-6:
433            #print 'h_C[%d]= %15.5e\n' % (i,h_C[i])
434            h_C[i] = 0.0
435            w_C[i] = z_C[i]
436            uh_C[i] = 0.0
437           
438    for name in [ 'stage', 'xmomentum']:
439        Q = domain.quantities[name]
440        if domain.order == 1:
441            Q.extrapolate_first_order()
442        elif domain.order == 2:
443            Q.extrapolate_second_order()
444        else:
445            raise 'Unknown order'
446
447    w_V  = domain.quantities['stage'].vertex_values                     
448    z_V  = domain.quantities['elevation'].vertex_values
449    h_V  = domain.quantities['height'].vertex_values
450    u_V  = domain.quantities['velocity'].vertex_values         
451    uh_V = domain.quantities['xmomentum'].vertex_values
452
453    h_V[:]    = w_V - z_V
454
455    for i in range(len(h_C)):
456        for j in range(2):
457            if h_V[i,j] < 0.0 :
458                #print 'h_V[%d,%d] = %f \n' % (i,j,h_V[i,j])                 
459                dh = h_V[i,j]
460                h_V[i,j] = 0.0
461                w_V[i,j] = z_V[i,j]
462                h_V[i,(j+1)%2] = h_V[i,(j+1)%2] + dh
463                w_V[i,(j+1)%2] = w_V[i,(j+1)%2] + dh
464                u_V[i,j] = 0.0
465            if h_V[i,j] < h0:
466                u_V[i,j]
467            else:
468                u_V[i,j] = uh_V[i,j]/(h_V[i,j] + h0/h_V[i,j])
469
470#---------------------------------------------------------------------------
471def distribute_to_vertices_and_edges_limit_a_d(domain):
472    """Distribution from centroids to vertices specific to the
473    shallow water wave equation.
474
475    In addition, all conserved quantities get distributed as per either a
476    constant (order==1) or a piecewise linear function (order==2).
477
478    Precondition:
479      All quantities defined at centroids and bed elevation defined at
480      vertices.
481
482    Postcondition
483      Conserved quantities defined at vertices
484
485    """
486
487    import sys
488    from Numeric import zeros, Float
489    from config import epsilon, h0
490       
491    N = domain.number_of_elements
492
493    #Shortcuts
494    Area = domain.quantities['area']
495    Discharge = domain.quantities['discharge']
496    Bed = domain.quantities['elevation']
497    Height = domain.quantities['height']
498    Stage = domain.quantities['stage']
499
500    #Arrays   
501    a_C   = Area.centroid_values   
502    d_C   = Discharge.centroid_values   
503    z_C   = Bed.centroid_values
504    h_C   = Height.centroid_values
505    w_C   = Stage.centroid_values
506       
507#work out stage
508    for i in range(N):
509        h_C[i] = w_C[i] - z_C[i]
510#make sure depth isn't zero       
511        if h_C[i] <= 1.0e-6:
512            #print 'h_C[%d]= %15.5e\n' % (i,h_C[i])
513            h_C[i] = 0.0
514            w_C[i] = z_C[i]
515            d_C[i] = 0.0
516#distribute stage,discharge
517
518    for name in [ 'area', 'discharge']:
519        Q = domain.quantities[name]
520        if domain.order == 1:
521            Q.extrapolate_first_order()
522        elif domain.order == 2:
523            Q.extrapolate_second_order()
524        else:
525            raise 'Unknown order'
526
527    a_V  = domain.quantities['area'].vertex_values                     
528    d_V  = domain.quantities['discharge'].vertex_values
529    h_V  = domain.quantities['height'].vertex_values
530    z_V  = domain.quantities['elevation'].vertex_values         
531    w_V = domain.quantities['stage'].vertex_values     
532#height at verticies
533    h_V[:]    = w_V - z_V
534
535##     for i in range(len(h_C)):
536##         for j in range(2):
537##             if h_V[i,j] < 0.0 :
538##                 #print 'h_V[%d,%d] = %f \n' % (i,j,h_V[i,j])                 
539##                 dh = h_V[i,j]
540##                 h_V[i,j] = 0.0
541##                 w_V[i,j] = z_V[i,j]
542##                 h_V[i,(j+1)%2] = h_V[i,(j+1)%2] + dh
543##                 w_V[i,(j+1)%2] = w_V[i,(j+1)%2] + dh
544##                 u_V[i,j] = 0.0
545##             if h_V[i,j] < h0:
546##                 u_V[i,j]
547##             else:
548##                 u_V[i,j] = uh_V[i,j]/(h_V[i,j] + h0/h_V[i,j]
549##                                        )
550               
551#--------------------------------------------------------
552#Boundaries - specific to the shallow_water_vel_domain
553#--------------------------------------------------------
554class Reflective_boundary(Boundary):
555    """Reflective boundary returns same conserved quantities as
556    those present in its neighbour volume but reflected.
557
558    This class is specific to the shallow water equation as it
559    works with the momentum quantities assumed to be the second
560    and third conserved quantities.
561    """
562
563    def __init__(self, domain = None):
564        Boundary.__init__(self)
565
566        if domain is None:
567            msg = 'Domain must be specified for reflective boundary'
568            raise msg
569
570        #Handy shorthands
571        self.normals  = domain.normals
572        self.stage    = domain.quantities['stage'].vertex_values
573        self.xmom     = domain.quantities['xmomentum'].vertex_values
574        self.bed      = domain.quantities['elevation'].vertex_values
575        self.height   = domain.quantities['height'].vertex_values
576        self.velocity = domain.quantities['velocity'].vertex_values
577
578        from Numeric import zeros, Float
579        #self.conserved_quantities = zeros(3, Float)
580        self.evolved_quantities = zeros(5, Float)
581
582    def __repr__(self):
583        return 'Reflective_boundary'
584
585
586    def evaluate(self, vol_id, edge_id):
587        """Reflective boundaries reverses the outward momentum
588        of the volume they serve.
589        """
590#Commenting out some quantities not currently keeping track of
591##         q = self.evolved_quantities
592##         q[0] = self.stage[vol_id, edge_id]
593##         q[1] = -self.xmom[vol_id, edge_id]
594##         q[2] = self.bed[vol_id, edge_id]
595##         q[3] = self.height[vol_id, edge_id]
596##         q[4] = -self.stage[stage_id, stage_id]
597
598        #print "In Reflective q ",q
599
600
601        return q
602
603class Dirichlet_boundary(Boundary):
604    """Dirichlet boundary returns constant values for the
605    conserved quantities
606    """
607
608
609    def __init__(self, evolved_quantities=None):
610        Boundary.__init__(self)
611
612        if evolved_quantities is None:
613            msg = 'Must specify one value for each evolved quantity'
614            raise msg
615
616        from Numeric import array, Float
617        self.evolved_quantities=array(evolved_quantities).astype(Float)
618
619    def __repr__(self):
620        return 'Dirichlet boundary (%s)' %self.evolved_quantities
621
622    def evaluate(self, vol_id=None, edge_id=None):
623        return self.evolved_quantities
624
625#--------------------------------------------------------
626#Boundaries for channel - specific to the channel domain
627#--------------------------------------------------------
628class Reflective_boundary(Boundary):
629    """Reflective boundary returns same conserved quantities as
630    those present in its neighbour volume but reflected.
631
632    This class is specific to the shallow water equation as it
633    works with the momentum quantities assumed to be the second
634    and third conserved quantities.
635    """
636
637    def __init__(self, domain = None):
638        Boundary.__init__(self)
639
640        if domain is None:
641            msg = 'Domain must be specified for reflective boundary'
642            raise msg
643
644        #Handy shorthands
645        self.normals  = domain.normals
646        self.area    = domain.quantities['area'].vertex_values
647        self.discharge     = domain.quantities['discharge'].vertex_values
648        self.bed      = domain.quantities['elevation'].vertex_values
649        self.height   = domain.quantities['height'].vertex_values
650        self.stage = domain.quantities['stage'].vertex_values
651
652        from Numeric import zeros, Float
653        #self.conserved_quantities = zeros(3, Float)
654        self.evolved_quantities = zeros(5, Float)
655
656    def __repr__(self):
657        return 'Reflective_boundary'
658
659
660    def evaluate(self, vol_id, edge_id):
661        """Reflective boundaries reverses the outward momentum
662        of the volume they serve.
663        """
664
665        q = self.evolved_quantities
666        q[0] = self.area[vol_id, edge_id]
667        q[1] = -self.discharge[vol_id, edge_id]
668        q[2] = self.bed[vol_id, edge_id]
669        q[3] = self.height[vol_id, edge_id]
670        q[4] = self.stage[vol_id, edge_id]
671
672        #print "In Reflective q ",q
673
674
675        return q
676
677class Dirichlet_boundary(Boundary):
678    """Dirichlet boundary returns constant values for the
679    conserved quantities
680    """
681
682
683    def __init__(self, evolved_quantities=None):
684        Boundary.__init__(self)
685
686        if evolved_quantities is None:
687            msg = 'Must specify one value for each evolved quantity'
688            raise msg
689
690        from Numeric import array, Float
691        self.evolved_quantities=array(evolved_quantities).astype(Float)
692
693    def __repr__(self):
694        return 'Dirichlet boundary (%s)' %self.evolved_quantities
695
696    def evaluate(self, vol_id=None, edge_id=None):
697        return self.evolved_quantities
698   
699#----------------------------
700#Standard forcing terms:
701#---------------------------
702def gravity(domain):
703    """Apply gravitational pull in the presence of bed slope
704    """
705
706    from util import gradient
707    from Numeric import zeros, Float, array, sum
708
709
710
711    Area     = domain.quantities['area']
712    Discharge      = domain.quantities['discharge']
713    Elevation = domain.quantities['elevation']
714    Height    = domain.quantities['height']
715    Stage     = domain.quantities['stage']
716
717    discharge_ud  = Discharge.explicit_update
718    #stage_ud = Stage.explicit_update
719
720
721    #h = Stage.vertex_values - Elevation.vertex_values
722    h = Height.vertex_values
723    b = Elevation.vertex_values
724    w = Stage.vertex_values
725
726    x = domain.get_vertex_coordinates()
727    g = domain.g
728
729    for k in range(domain.number_of_elements):
730        avg_h = 0.5*(h[k,0] + h[k,1])
731
732        #Compute bed slope
733        x0, x1 = x[k,:]
734        b0, b1 = b[k,:]
735        bx = gradient(x0, x1, b0, b1)
736       
737        #Update momentum (explicit update is reset to source values)
738        discharge_ud[k] += -g*bx*avg_h
739        #stage_ud[k] = 0.0
740 
741 
742def manning_friction(domain):
743    """Apply (Manning) friction to water momentum
744    """
745
746    from math import sqrt
747
748    w = domain.quantities['stage'].centroid_values
749    z = domain.quantities['elevation'].centroid_values
750    h = w-z
751
752    uh = domain.quantities['xmomentum'].centroid_values
753    #vh = domain.quantities['ymomentum'].centroid_values
754    eta = domain.quantities['friction'].centroid_values
755
756    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
757    #ymom_update = domain.quantities['ymomentum'].semi_implicit_update
758
759    N = domain.number_of_elements
760    eps = domain.minimum_allowed_height
761    g = domain.g
762
763    for k in range(N):
764        if eta[k] >= eps:
765            if h[k] >= eps:
766                #S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
767                S = -g * eta[k]**2 * uh[k]
768                S /= h[k]**(7.0/3)
769
770                #Update momentum
771                xmom_update[k] += S*uh[k]
772                #ymom_update[k] += S*vh[k]
773
774def linear_friction(domain):
775    """Apply linear friction to water momentum
776
777    Assumes quantity: 'linear_friction' to be present
778    """
779
780    from math import sqrt
781
782    w = domain.quantities['stage'].centroid_values
783    z = domain.quantities['elevation'].centroid_values
784    h = w-z
785
786    uh = domain.quantities['xmomentum'].centroid_values
787    tau = domain.quantities['linear_friction'].centroid_values
788
789    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
790
791    N = domain.number_of_elements
792    eps = domain.minimum_allowed_height
793
794    for k in range(N):
795        if tau[k] >= eps:
796            if h[k] >= eps:
797                S = -tau[k]/h[k]
798
799                #Update momentum
800                xmom_update[k] += S*uh[k]
801
802
803
804def check_forcefield(f):
805    """Check that f is either
806    1: a callable object f(t,x,y), where x and y are vectors
807       and that it returns an array or a list of same length
808       as x and y
809    2: a scalar
810    """
811
812    from Numeric import ones, Float, array
813
814
815    if callable(f):
816        #N = 3
817        N = 2
818        #x = ones(3, Float)
819        #y = ones(3, Float)
820        x = ones(2, Float)
821        #y = ones(2, Float)
822       
823        try:
824            #q = f(1.0, x=x, y=y)
825            q = f(1.0, x=x)
826        except Exception, e:
827            msg = 'Function %s could not be executed:\n%s' %(f, e)
828            #FIXME: Reconsider this semantics
829            raise msg
830
831        try:
832            q = array(q).astype(Float)
833        except:
834            msg = 'Return value from vector function %s could ' %f
835            msg += 'not be converted into a Numeric array of floats.\n'
836            msg += 'Specified function should return either list or array.'
837            raise msg
838
839        #Is this really what we want?
840        msg = 'Return vector from function %s ' %f
841        msg += 'must have same lenght as input vectors'
842        assert len(q) == N, msg
843
844    else:
845        try:
846            f = float(f)
847        except:
848            msg = 'Force field %s must be either a scalar' %f
849            msg += ' or a vector function'
850            raise msg
851    return f
852
Note: See TracBrowser for help on using the repository browser.