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

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