source: anuga_work/development/anuga_1d/shallow_water_domain_new.py @ 5741

Last change on this file since 5741 was 5741, checked in by steve, 16 years ago

Updated timestepping

File size: 28.7 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        other_quantities = ['elevation', 'height', 'velocity']
56        Generic_Domain.__init__(self, coordinates, boundary,
57                                conserved_quantities, other_quantities,
58                                tagged_elements)
59       
60        from config import minimum_allowed_height, g
61        self.minimum_allowed_height = minimum_allowed_height
62        self.g = g
63
64        #forcing terms not included in 1d domain ? Why?
65        self.forcing_terms.append(gravity)
66       
67        #Realtime visualisation
68        self.visualiser = None
69        self.visualise  = False
70        self.visualise_color_stage = False
71        self.visualise_stage_range = 1.0
72        self.visualise_timer = True
73        self.visualise_range_z = None
74       
75        #Stored output
76        self.store = True
77        self.format = 'sww'
78        self.smooth = True
79
80        #Evolve parametrs
81        self.cfl = 1.0
82       
83        #Reduction operation for get_vertex_values
84        from util import mean
85        self.reduction = mean
86        #self.reduction = min  #Looks better near steep slopes
87
88        self.quantities_to_be_stored = ['stage', 'xmomentum', 'elevation']
89
90        self.__doc__ = 'shallow_water_domain_new'
91
92
93    def set_quantities_to_be_stored(self, q):
94        """Specify which quantities will be stored in the sww file.
95
96        q must be either:
97          - the name of a quantity
98          - a list of quantity names
99          - None
100
101        In the two first cases, the named quantities will be stored at each
102        yieldstep
103        (This is in addition to the quantities elevation and friction) 
104        If q is None, storage will be switched off altogether.
105        """
106
107
108        if q is None:
109            self.quantities_to_be_stored = []           
110            self.store = False
111            return
112
113        if isinstance(q, basestring):
114            q = [q] # Turn argument into a list
115
116        #Check correcness   
117        for quantity_name in q:
118            msg = 'Quantity %s is not a valid conserved quantity' %quantity_name
119            assert quantity_name in self.conserved_quantities, msg
120       
121        self.quantities_to_be_stored = q
122       
123
124    def initialise_visualiser(self,scale_z=1.0,rect=None):
125        #Realtime visualisation
126        if self.visualiser is None:
127            from realtime_visualisation_new import Visualiser
128            self.visualiser = Visualiser(self,scale_z,rect)
129            self.visualiser.setup['elevation']=True
130            self.visualiser.updating['stage']=True
131        self.visualise = True
132        if self.visualise_color_stage == True:
133            self.visualiser.coloring['stage'] = True
134            self.visualiser.qcolor['stage'] = (0.0, 0.0, 0.8)
135        print 'initialise visualiser'
136        print self.visualiser.setup
137        print self.visualiser.updating
138
139    def check_integrity(self):
140        Generic_Domain.check_integrity(self)
141        #Check that we are solving the shallow water wave equation
142
143        msg = 'First conserved quantity must be "stage"'
144        assert self.conserved_quantities[0] == 'stage', msg
145        msg = 'Second conserved quantity must be "xmomentum"'
146        assert self.conserved_quantities[1] == 'xmomentum', msg
147
148    def extrapolate_second_order_sw(self):
149        #Call correct module function
150        #(either from this module or C-extension)
151        extrapolate_second_order_sw(self)
152
153    def compute_fluxes(self):
154        #Call correct module function
155        #(either from this module or C-extension)
156        compute_fluxes_C_short(self)    #compute_fluxes_C_wellbalanced(self)     #compute_fluxes(self)  #compute_fluxes_C_long(self)
157
158    def compute_timestep(self):
159        #Call correct module function
160        compute_timestep(self)
161
162    def distribute_to_vertices_and_edges(self):
163        #Call correct module function
164        #(either from this module or C-extension)
165        distribute_to_vertices_and_edges(self)
166       
167    def evolve(self, yieldstep = None,
168               finaltime = None,
169               duration = None,
170               skip_initial_step = False):
171        """Specialisation of basic evolve method from parent class
172        """
173
174        self.distribute_to_vertices_and_edges()
175       
176        for t in Generic_Domain.evolve(self, yieldstep, finaltime, duration,
177                                       skip_initial_step):
178            #Pass control on to outer loop for more specific actions
179            yield(t)
180
181    def initialise_storage(self):
182        """Create and initialise self.writer object for storing data.
183        Also, save x and bed elevation
184        """
185
186        import data_manager
187
188        #Initialise writer
189        self.writer = data_manager.get_dataobject(self, mode = 'w')
190
191        #Store vertices and connectivity
192        self.writer.store_connectivity()
193
194
195    def store_timestep(self, name):
196        """Store named quantity and time.
197
198        Precondition:
199           self.write has been initialised
200        """
201        self.writer.store_timestep(name)
202
203
204#=============== End of Shallow Water Domain ===============================
205def compute_timestep(domain):
206    import sys
207    from Numeric import zeros, Float
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   
216    stage = Stage.vertex_values
217    xmom =  Xmom.vertex_values
218    bed =   Bed.vertex_values
219
220    stage_bdry = Stage.boundary_values
221    xmom_bdry =  Xmom.boundary_values
222
223    flux = zeros(2, Float) #Work array for summing up fluxes
224    ql = zeros(2, Float)
225    qr = zeros(2, Float)
226
227    #Loop
228    timestep = float(sys.maxint)
229    enter = True
230    for k in range(N):
231
232        flux[:] = 0.  #Reset work array
233        for i in range(2):
234            #Quantities inside volume facing neighbour i
235            ql = [stage[k, i], xmom[k, i]]
236            zl = bed[k, i]
237
238            #Quantities at neighbour on nearest face
239            n = domain.neighbours[k,i]
240            if n < 0:
241                m = -n-1 #Convert negative flag to index
242                qr[0] = stage_bdry[m]
243                qr[1] = xmom_bdry[m]
244                zr = zl #Extend bed elevation to boundary
245            else:
246                #m = domain.neighbour_edges[k,i]
247                m = domain.neighbour_vertices[k,i]
248                qr[0] = stage[n, m]
249                qr[1] =  xmom[n, m]
250                zr = bed[n, m]
251
252
253            #Outward pointing normal vector
254            normal = domain.normals[k, i]
255
256            #if domain.split == False:
257            #    edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
258            #elif domain.split == True:
259            #    edgeflux, max_speed = flux_function_split(normal, ql, qr, zl, zr)
260            #Update optimal_timestep
261            #try:
262            #    timestep = min(timestep, domain.cfl*0.5*domain.areas[k]/max_speed)
263            #except ZeroDivisionError:
264            #    pass
265
266    domain.timestep = timestep
267#
268def flux_function(normal, ql, qr, zl, zr):
269    """Compute fluxes between volumes for the shallow water wave equation
270    cast in terms of w = h+z using the 'central scheme' as described in
271
272    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
273    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
274    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
275
276    The implemented formula is given in equation (3.15) on page 714
277
278    Conserved quantities w, uh, are stored as elements 0 and 1
279    in the numerical vectors ql an qr.
280
281    Bed elevations zl and zr.
282    """
283
284    from config import g, epsilon
285    from math import sqrt
286    from Numeric import array
287
288    #print 'ql',ql
289
290    #Align momentums with x-axis
291    q_left = ql
292    q_left[1] = q_left[1]*normal
293    q_right = qr
294    q_right[1] = q_right[1]*normal
295
296    z = 0.5*(zl+zr) #Take average of field values
297
298    w_left  = q_left[0]   #w=h+z
299    h_left  = w_left-z
300    uh_left = q_left[1]
301
302    if h_left < epsilon:
303        u_left = 0.0  #Could have been negative
304        h_left = 0.0
305    else:
306        u_left  = uh_left/h_left
307
308
309    w_right  = q_right[0]  #w=h+z
310    h_right  = w_right-z
311    uh_right = q_right[1]
312
313
314    if h_right < epsilon:
315        u_right = 0.0  #Could have been negative
316        h_right = 0.0
317    else:
318        u_right  = uh_right/h_right
319
320    soundspeed_left  = sqrt(g*h_left)
321    soundspeed_right = sqrt(g*h_right)
322
323    #Maximal wave speed
324    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
325   
326    #Minimal wave speed
327    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
328   
329    #Flux computation
330    flux_left  = array([u_left*h_left,
331                        u_left*uh_left + 0.5*g*h_left*h_left])
332    flux_right = array([u_right*h_right,
333                        u_right*uh_right + 0.5*g*h_right*h_right])
334
335    denom = s_max-s_min
336    if denom == 0.0:
337        edgeflux = array([0.0, 0.0])
338        max_speed = 0.0
339    else:
340        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
341        edgeflux += s_max*s_min*(q_right-q_left)/denom
342       
343        edgeflux[1] = edgeflux[1]*normal
344
345        max_speed = max(abs(s_max), abs(s_min))
346
347    return edgeflux, max_speed
348
349#The compute_fluxes is involving the flux_function above.
350def compute_fluxes(domain):
351    """Compute all fluxes and the timestep suitable for all volumes
352    in domain.
353
354    Compute total flux for each conserved quantity using "flux_function"
355
356    Fluxes across each edge are scaled by edgelengths and summed up
357    Resulting flux is then scaled by area and stored in
358    explicit_update for each of the three conserved quantities
359    stage, xmomentum and ymomentum
360
361    The maximal allowable speed computed by the flux_function for each volume
362    is converted to a timestep that must not be exceeded. The minimum of
363    those is computed as the next overall timestep.
364
365    Post conditions:
366      domain.explicit_update is reset to computed flux values
367      domain.timestep is set to the largest step satisfying all volumes.
368    """
369
370    import sys
371    from Numeric import zeros, Float
372
373    N = domain.number_of_elements
374
375    #Shortcuts
376    Stage = domain.quantities['stage']
377    Xmom = domain.quantities['xmomentum']
378#    Ymom = domain.quantities['ymomentum']
379    Bed = domain.quantities['elevation']
380
381    #Arrays
382    stage = Stage.vertex_values
383    xmom =  Xmom.vertex_values
384    bed =   Bed.vertex_values
385   
386    stage_bdry = Stage.boundary_values
387    xmom_bdry =  Xmom.boundary_values
388   
389    flux = zeros(2, Float) #Work array for summing up fluxes
390    ql = zeros(2, Float)
391    qr = zeros(2, Float)
392
393    #Loop
394    timestep = float(sys.maxint)
395    enter = True
396    for k in range(N):
397
398        flux[:] = 0.  #Reset work array
399        #for i in range(3):
400        for i in range(2):
401            #Quantities inside volume facing neighbour i
402            ql = [stage[k, i], xmom[k, i]]
403            zl = bed[k, i]
404
405            #Quantities at neighbour on nearest face
406            n = domain.neighbours[k,i]
407            if n < 0:
408                m = -n-1 #Convert negative flag to index
409                qr[0] = stage_bdry[m]
410                qr[1] = xmom_bdry[m]
411                zr = zl #Extend bed elevation to boundary
412            else:
413                #m = domain.neighbour_edges[k,i]
414                m = domain.neighbour_vertices[k,i]
415                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
416                qr[0] = stage[n, m]
417                qr[1] = xmom[n, m]
418                zr = bed[n, m]
419
420            #Outward pointing normal vector
421            normal = domain.normals[k, i]
422       
423            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
424
425            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
426            # flux = edgefluxleft - edgefluxright
427            flux -= edgeflux #* domain.edgelengths[k,i]
428            #Update optimal_timestep
429            try:
430                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
431                timestep = min(timestep, domain.cfl*0.5*domain.areas[k]/max_speed)
432            except ZeroDivisionError:
433                pass
434
435        #Normalise by area and store for when all conserved
436        #quantities get updated
437        flux /= domain.areas[k]
438
439        Stage.explicit_update[k] = flux[0]
440        Xmom.explicit_update[k] =  flux[1]
441       
442
443    domain.flux_timestep = timestep
444    #print domain.quantities['stage'].centroid_values
445
446
447# Compute flux definition
448def compute_fluxes_C_long(domain):
449        from Numeric import zeros, Float
450        import sys
451
452       
453        timestep = float(sys.maxint)
454        epsilon = domain.epsilon
455        g = domain.g
456        neighbours = domain.neighbours
457        neighbour_vertices = domain.neighbour_vertices
458        normals = domain.normals
459        areas = domain.areas
460        stage_edge_values = domain.quantities['stage'].vertex_values
461        xmom_edge_values = domain.quantities['xmomentum'].vertex_values
462        bed_edge_values = domain.quantities['elevation'].vertex_values
463        stage_boundary_values = domain.quantities['stage'].boundary_values
464        xmom_boundary_values = domain.quantities['xmomentum'].boundary_values
465        stage_explicit_update = domain.quantities['stage'].explicit_update
466        xmom_explicit_update = domain.quantities['xmomentum'].explicit_update
467        number_of_elements = len(stage_edge_values)
468        max_speed_array = domain.max_speed_array
469               
470        from comp_flux_ext import compute_fluxes_ext
471       
472        domain.flux_timestep = compute_fluxes_ext(timestep,
473                                  epsilon,
474                                  g,
475                                  neighbours,
476                                  neighbour_vertices,
477                                  normals,
478                                  areas, 
479                                  stage_edge_values,
480                                  xmom_edge_values,
481                                  bed_edge_values,
482                                  stage_boundary_values,
483                                  xmom_boundary_values,
484                                  stage_explicit_update,
485                                  xmom_explicit_update,
486                                  number_of_elements,
487                                  max_speed_array)
488
489
490# Compute flux definition
491def compute_fluxes_C_short(domain):
492        from Numeric import zeros, Float
493        import sys
494
495       
496        timestep = float(sys.maxint)
497
498        stage = domain.quantities['stage']
499        xmom = domain.quantities['xmomentum']
500        bed = domain.quantities['elevation']
501
502
503        from comp_flux_ext import compute_fluxes_ext_short
504
505        domain.flux_timestep = compute_fluxes_ext_short(timestep,domain,stage,xmom,bed)
506
507
508def compute_fluxes_C_wellbalanced(domain):
509        from Numeric import zeros, Float
510        import sys
511
512       
513        timestep = float(sys.maxint)
514        epsilon = domain.epsilon
515        g = domain.g
516        neighbours = domain.neighbours
517        neighbour_vertices = domain.neighbour_vertices
518        normals = domain.normals
519        areas = domain.areas
520        stage_edge_values = domain.quantities['stage'].vertex_values
521        xmom_edge_values = domain.quantities['xmomentum'].vertex_values
522        bed_edge_values = domain.quantities['elevation'].vertex_values
523        stage_boundary_values = domain.quantities['stage'].boundary_values
524        xmom_boundary_values = domain.quantities['xmomentum'].boundary_values
525        stage_explicit_update = domain.quantities['stage'].explicit_update
526        xmom_explicit_update = domain.quantities['xmomentum'].explicit_update
527        number_of_elements = len(stage_edge_values)
528        max_speed_array = domain.max_speed_array
529               
530        from comp_flux_ext_wellbalanced import compute_fluxes_ext_wellbalanced  #from comp_flux_ext import compute_fluxes_ext
531       
532        domain.flux_timestep = compute_fluxes_ext_wellbalanced(timestep,
533                                  epsilon,
534                                  g,
535                                  neighbours,
536                                  neighbour_vertices,
537                                  normals,
538                                  areas, 
539                                  stage_edge_values,
540                                  xmom_edge_values,
541                                  bed_edge_values,
542                                  stage_boundary_values,
543                                  xmom_boundary_values,
544                                  stage_explicit_update,
545                                  xmom_explicit_update,
546                                  number_of_elements,
547                                  max_speed_array)
548
549# ###################################
550
551
552
553
554
555
556# Module functions for gradient limiting (distribute_to_vertices_and_edges)
557
558def distribute_to_vertices_and_edges(domain):
559    """Distribution from centroids to vertices specific to the
560    shallow water wave
561    equation.
562
563    It will ensure that h (w-z) is always non-negative even in the
564    presence of steep bed-slopes by taking a weighted average between shallow
565    and deep cases.
566
567    In addition, all conserved quantities get distributed as per either a
568    constant (order==1) or a piecewise linear function (order==2).
569
570    FIXME: more explanation about removal of artificial variability etc
571
572    Precondition:
573      All quantities defined at centroids and bed elevation defined at
574      vertices.
575
576    Postcondition
577      Conserved quantities defined at vertices
578
579    """
580
581    #from config import optimised_gradient_limiter
582
583    #Remove very thin layers of water
584    protect_against_infinitesimal_and_negative_heights(domain)
585   
586
587    for name in domain.conserved_quantities:
588        Q = domain.quantities[name]
589        if domain.order == 1:
590            Q.extrapolate_first_order()
591        elif domain.order == 2:
592            #print "add extrapolate second order to shallow water"
593            #if name != 'height':
594            Q.extrapolate_second_order()
595            #Q.limit()
596        else:
597            raise 'Unknown order'
598
599   
600def protect_against_infinitesimal_and_negative_heights(domain):
601    """Protect against infinitesimal heights and associated high velocities
602    """
603
604    #Shortcuts
605    wc = domain.quantities['stage'].centroid_values
606    zc = domain.quantities['elevation'].centroid_values
607    xmomc = domain.quantities['xmomentum'].centroid_values
608#    ymomc = domain.quantities['ymomentum'].centroid_values
609    hc = wc - zc  #Water depths at centroids
610
611    zv = domain.quantities['elevation'].vertex_values
612    wv = domain.quantities['stage'].vertex_values
613    hv = wv-zv
614    xmomv = domain.quantities['xmomentum'].vertex_values
615    #remove the above two lines and corresponding code below
616
617    #Update
618    for k in range(domain.number_of_elements):
619
620        if hc[k] < domain.minimum_allowed_height:
621            #Control stage
622            if hc[k] < domain.epsilon:
623                wc[k] = zc[k] # Contain 'lost mass' error
624                wv[k,0] = zv[k,0]
625                wv[k,1] = zv[k,1]
626               
627            xmomc[k] = 0.0
628           
629       
630   
631def h_limiter(domain):
632    """Limit slopes for each volume to eliminate artificial variance
633    introduced by e.g. second order extrapolator
634
635    limit on h = w-z
636
637    This limiter depends on two quantities (w,z) so it resides within
638    this module rather than within quantity.py
639    """
640
641    from Numeric import zeros, Float
642
643    N = domain.number_of_elements
644    beta_h = domain.beta_h
645
646    #Shortcuts
647    wc = domain.quantities['stage'].centroid_values
648    zc = domain.quantities['elevation'].centroid_values
649    hc = wc - zc
650
651    wv = domain.quantities['stage'].vertex_values
652    zv = domain.quantities['elevation'].vertex_values
653    hv = wv-zv
654
655    hvbar = zeros(hv.shape, Float) #h-limited values
656
657    #Find min and max of this and neighbour's centroid values
658    hmax = zeros(hc.shape, Float)
659    hmin = zeros(hc.shape, Float)
660
661    for k in range(N):
662        hmax[k] = hmin[k] = hc[k]
663        #for i in range(3):
664        for i in range(2):   
665            n = domain.neighbours[k,i]
666            if n >= 0:
667                hn = hc[n] #Neighbour's centroid value
668
669                hmin[k] = min(hmin[k], hn)
670                hmax[k] = max(hmax[k], hn)
671
672
673    #Diffences between centroids and maxima/minima
674    dhmax = hmax - hc
675    dhmin = hmin - hc
676
677    #Deltas between vertex and centroid values
678    dh = zeros(hv.shape, Float)
679    #for i in range(3):
680    for i in range(2):
681        dh[:,i] = hv[:,i] - hc
682
683    #Phi limiter
684    for k in range(N):
685
686        #Find the gradient limiter (phi) across vertices
687        phi = 1.0
688        #for i in range(3):
689        for i in range(2):
690            r = 1.0
691            if (dh[k,i] > 0): r = dhmax[k]/dh[k,i]
692            if (dh[k,i] < 0): r = dhmin[k]/dh[k,i]
693
694            phi = min( min(r*beta_h, 1), phi )
695
696        #Then update using phi limiter
697        #for i in range(3):
698        for i in range(2):
699            hvbar[k,i] = hc[k] + phi*dh[k,i]
700
701    return hvbar
702
703def balance_deep_and_shallow(domain):
704    """Compute linear combination between stage as computed by
705    gradient-limiters limiting using w, and stage computed by
706    gradient-limiters limiting using h (h-limiter).
707    The former takes precedence when heights are large compared to the
708    bed slope while the latter takes precedence when heights are
709    relatively small.  Anything in between is computed as a balanced
710    linear combination in order to avoid numerical disturbances which
711    would otherwise appear as a result of hard switching between
712    modes.
713
714    The h-limiter is always applied irrespective of the order.
715    """
716
717    #Shortcuts
718    wc = domain.quantities['stage'].centroid_values
719    zc = domain.quantities['elevation'].centroid_values
720    hc = wc - zc
721
722    wv = domain.quantities['stage'].vertex_values
723    zv = domain.quantities['elevation'].vertex_values
724    hv = wv-zv
725
726    #Limit h
727    hvbar = h_limiter(domain)
728
729    for k in range(domain.number_of_elements):
730       
731        dz = max(abs(zv[k,0]-zc[k]),
732                 abs(zv[k,1]-zc[k]))#,
733#                 abs(zv[k,2]-zc[k]))
734
735
736        hmin = min( hv[k,:] )
737
738        if dz > 0.0:
739            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
740        else:
741            #Flat bed
742            alpha = 1.0
743
744        alpha  = 0.0
745       
746        if alpha < 1:
747
748            #for i in range(3):
749            for i in range(2):
750                wv[k,i] = zv[k,i] + (1.0-alpha)*hvbar[k,i] + alpha*hv[k,i]
751
752            #Momentums at centroids
753            xmomc = domain.quantities['xmomentum'].centroid_values
754#            ymomc = domain.quantities['ymomentum'].centroid_values
755
756            #Momentums at vertices
757            xmomv = domain.quantities['xmomentum'].vertex_values
758#           ymomv = domain.quantities['ymomentum'].vertex_values
759
760            # Update momentum as a linear combination of
761            # xmomc and ymomc (shallow) and momentum
762            # from extrapolator xmomv and ymomv (deep).
763            xmomv[k,:] = (1.0-alpha)*xmomc[k] + alpha*xmomv[k,:]
764#            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
765
766
767# ##############################################
768#Boundaries - specific to the shallow water wave equation
769class Reflective_boundary(Boundary):
770    """Reflective boundary returns same conserved quantities as
771    those present in its neighbour volume but reflected.
772
773    This class is specific to the shallow water equation as it
774    works with the momentum quantities assumed to be the second
775    and third conserved quantities.
776    """
777
778    def __init__(self, domain = None):
779        Boundary.__init__(self)
780
781        if domain is None:
782            msg = 'Domain must be specified for reflective boundary'
783            raise msg
784
785        #Handy shorthands
786        self.normals = domain.normals
787        self.stage   = domain.quantities['stage'].vertex_values
788        self.xmom    = domain.quantities['xmomentum'].vertex_values
789
790        from Numeric import zeros, Float
791        #self.conserved_quantities = zeros(3, Float)
792        self.conserved_quantities = zeros(2, Float)
793
794    def __repr__(self):
795        return 'Reflective_boundary'
796
797
798    def evaluate(self, vol_id, edge_id):
799        """Reflective boundaries reverses the outward momentum
800        of the volume they serve.
801        """
802
803        q = self.conserved_quantities
804        q[0] = self.stage[vol_id, edge_id]
805        q[1] = self.xmom[vol_id, edge_id]
806        normal = self.normals[vol_id,edge_id]
807
808        r = q
809        r[1] = normal*r[1]
810        r[1] = -r[1]
811        r[1] = normal*r[1]
812        q = r
813        #For start interval there is no outward momentum so do not need to
814        #reverse direction in this case
815
816        return q
817
818class Dirichlet_boundary(Boundary):
819    """Dirichlet boundary returns constant values for the
820    conserved quantities
821    """
822
823
824    def __init__(self, conserved_quantities=None):
825        Boundary.__init__(self)
826
827        if conserved_quantities is None:
828            msg = 'Must specify one value for each conserved quantity'
829            raise msg
830
831        from Numeric import array, Float
832        self.conserved_quantities=array(conserved_quantities).astype(Float)
833
834    def __repr__(self):
835        return 'Dirichlet boundary (%s)' %self.conserved_quantities
836
837    def evaluate(self, vol_id=None, edge_id=None):
838        return self.conserved_quantities
839
840
841# ########################
842#Standard forcing terms:
843def gravity(domain):
844    """Apply gravitational pull in the presence of bed slope
845    """
846
847    from util import gradient
848    from Numeric import zeros, Float, array, sum
849
850    xmom = domain.quantities['xmomentum'].explicit_update
851    stage = domain.quantities['stage'].explicit_update
852#    ymom = domain.quantities['ymomentum'].explicit_update
853
854    Stage = domain.quantities['stage']
855    Elevation = domain.quantities['elevation']
856    #h = Stage.edge_values - Elevation.edge_values
857    h = Stage.vertex_values - Elevation.vertex_values
858    b = Elevation.vertex_values
859    w = Stage.vertex_values
860
861    x = domain.get_vertex_coordinates()
862    g = domain.g
863
864    for k in range(domain.number_of_elements):
865#        avg_h = sum( h[k,:] )/3
866        avg_h = sum( h[k,:] )/2
867
868        #Compute bed slope
869        #x0, y0, x1, y1, x2, y2 = x[k,:]
870        x0, x1 = x[k,:]
871        #z0, z1, z2 = v[k,:]
872        b0, b1 = b[k,:]
873
874        w0, w1 = w[k,:]
875        wx = gradient(x0, x1, w0, w1)
876
877        #zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
878        bx = gradient(x0, x1, b0, b1)
879       
880        #Update momentum (explicit update is reset to source values)
881        xmom[k] += -g*bx*avg_h
882        #xmom[k] = -g*bx*avg_h
883        #stage[k] = 0.0
884 
885 
886
887def check_forcefield(f):
888    """Check that f is either
889    1: a callable object f(t,x,y), where x and y are vectors
890       and that it returns an array or a list of same length
891       as x and y
892    2: a scalar
893    """
894
895    from Numeric import ones, Float, array
896
897
898    if callable(f):
899        #N = 3
900        N = 2
901        #x = ones(3, Float)
902        #y = ones(3, Float)
903        x = ones(2, Float)
904        #y = ones(2, Float)
905       
906        try:
907            #q = f(1.0, x=x, y=y)
908            q = f(1.0, x=x)
909        except Exception, e:
910            msg = 'Function %s could not be executed:\n%s' %(f, e)
911            #FIXME: Reconsider this semantics
912            raise msg
913
914        try:
915            q = array(q).astype(Float)
916        except:
917            msg = 'Return value from vector function %s could ' %f
918            msg += 'not be converted into a Numeric array of floats.\n'
919            msg += 'Specified function should return either list or array.'
920            raise msg
921
922        #Is this really what we want?
923        msg = 'Return vector from function %s ' %f
924        msg += 'must have same lenght as input vectors'
925        assert len(q) == N, msg
926
927    else:
928        try:
929            f = float(f)
930        except:
931            msg = 'Force field %s must be either a scalar' %f
932            msg += ' or a vector function'
933            raise msg
934    return f
Note: See TracBrowser for help on using the repository browser.