source: anuga_work/development/anuga_1d/shallow_water_domain.py @ 5562

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

Added unit test files for quantity and shallow water

File size: 37.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#from domain import *
46#from domain_order2 import *
47from domain import *
48Generic_Domain = Domain #Rename
49
50#Shallow water domain
51class Domain(Generic_Domain):
52
53    def __init__(self, coordinates, boundary = None, tagged_elements = None):
54
55        conserved_quantities = ['stage', 'xmomentum']
56        other_quantities = ['elevation', 'friction']
57        Generic_Domain.__init__(self, coordinates, boundary,
58                                conserved_quantities, other_quantities,
59                                tagged_elements)
60       
61        from config import minimum_allowed_height, g
62        self.minimum_allowed_height = minimum_allowed_height
63        self.g = g
64
65        #forcing terms not included in 1d domain ?WHy?
66        self.forcing_terms.append(gravity)
67        #self.forcing_terms.append(manning_friction)
68        #print "\nI have Removed forcing terms line 64 1dsw"
69
70        #Realtime visualisation
71        self.visualiser = None
72        self.visualise  = False
73        self.visualise_color_stage = False
74        self.visualise_stage_range = 1.0
75        self.visualise_timer = True
76        self.visualise_range_z = None
77       
78        #Stored output
79        self.store = True
80        self.format = 'sww'
81        self.smooth = True
82
83        #Evolve parametrs
84        self.cfl = 1.0
85       
86        #Reduction operation for get_vertex_values
87        from util import mean
88        self.reduction = mean
89        #self.reduction = min  #Looks better near steep slopes
90
91        self.quantities_to_be_stored = ['stage','xmomentum']
92
93
94    def set_quantities_to_be_stored(self, q):
95        """Specify which quantities will be stored in the sww file.
96
97        q must be either:
98          - the name of a quantity
99          - a list of quantity names
100          - None
101
102        In the two first cases, the named quantities will be stored at each
103        yieldstep
104        (This is in addition to the quantities elevation and friction) 
105        If q is None, storage will be switched off altogether.
106        """
107
108
109        if q is None:
110            self.quantities_to_be_stored = []           
111            self.store = False
112            return
113
114        if isinstance(q, basestring):
115            q = [q] # Turn argument into a list
116
117        #Check correcness   
118        for quantity_name in q:
119            msg = 'Quantity %s is not a valid conserved quantity' %quantity_name
120            assert quantity_name in self.conserved_quantities, msg
121       
122        self.quantities_to_be_stored = q
123       
124
125    def initialise_visualiser(self,scale_z=1.0,rect=None):
126        #Realtime visualisation
127        if self.visualiser is None:
128            from realtime_visualisation_new import Visualiser
129            self.visualiser = Visualiser(self,scale_z,rect)
130            self.visualiser.setup['elevation']=True
131            self.visualiser.updating['stage']=True
132        self.visualise = True
133        if self.visualise_color_stage == True:
134            self.visualiser.coloring['stage'] = True
135            self.visualiser.qcolor['stage'] = (0.0, 0.0, 0.8)
136        print 'initialise visualiser'
137        print self.visualiser.setup
138        print self.visualiser.updating
139
140    def check_integrity(self):
141        Generic_Domain.check_integrity(self)
142        #Check that we are solving the shallow water wave equation
143
144        msg = 'First conserved quantity must be "stage"'
145        assert self.conserved_quantities[0] == 'stage', msg
146        msg = 'Second conserved quantity must be "xmomentum"'
147        assert self.conserved_quantities[1] == 'xmomentum', msg
148
149    def extrapolate_second_order_sw(self):
150        #Call correct module function
151        #(either from this module or C-extension)
152        extrapolate_second_order_sw(self)
153
154    def compute_fluxes(self):
155        #Call correct module function
156        #(either from this module or C-extension)
157        compute_fluxes(self)
158
159    def compute_timestep(self):
160        #Call correct module function
161        compute_timestep(self)
162
163    def distribute_to_vertices_and_edges(self):
164        #Call correct module function
165        #(either from this module or C-extension)
166        distribute_to_vertices_and_edges(self)
167       
168    def evolve(self, yieldstep = None, finaltime = None,
169               skip_initial_step = False):
170        """Specialisation of basic evolve method from parent class
171        """
172
173        #Call check integrity here rather than from user scripts
174        #self.check_integrity()
175
176        #msg = 'Parameter beta_h must be in the interval [0, 1)'
177        #assert 0 <= self.beta_h < 1.0, msg
178        #msg = 'Parameter beta_w must be in the interval [0, 1)'
179        #assert 0 <= self.beta_w < 1.0, msg
180
181
182        #Initial update of vertex and edge values before any storage
183        #and or visualisation
184       
185        self.distribute_to_vertices_and_edges()
186       
187        #Initialise real time viz if requested
188        #if self.visualise is True and self.time == 0.0:
189        #    if self.visualiser is None:
190        #        self.initialise_visualiser()
191        #
192        #    self.visualiser.update_timer()
193        #    self.visualiser.setup_all()
194
195        #Store model data, e.g. for visualisation
196        #if self.store is True and self.time == 0.0:
197        #    self.initialise_storage()
198        #    #print 'Storing results in ' + self.writer.filename
199        #else:
200        #    pass
201        #    #print 'Results will not be stored.'
202        #    #print 'To store results set domain.store = True'
203        #    #FIXME: Diagnostic output should be controlled by
204        #    # a 'verbose' flag living in domain (or in a parent class)
205
206        #Call basic machinery from parent class
207        for t in Generic_Domain.evolve(self, yieldstep, finaltime,
208                                       skip_initial_step):
209            #Real time viz
210        #    if self.visualise is True:
211        #        self.visualiser.update_all()
212        #        self.visualiser.update_timer()
213
214
215            #Store model data, e.g. for subsequent visualisation
216        #    if self.store is True:
217        #        self.store_timestep(self.quantities_to_be_stored)
218
219            #FIXME: Could maybe be taken from specified list
220            #of 'store every step' quantities
221
222            #Pass control on to outer loop for more specific actions
223            yield(t)
224
225    def initialise_storage(self):
226        """Create and initialise self.writer object for storing data.
227        Also, save x and bed elevation
228        """
229
230        import data_manager
231
232        #Initialise writer
233        self.writer = data_manager.get_dataobject(self, mode = 'w')
234
235        #Store vertices and connectivity
236        self.writer.store_connectivity()
237
238
239    def store_timestep(self, name):
240        """Store named quantity and time.
241
242        Precondition:
243           self.write has been initialised
244        """
245        self.writer.store_timestep(name)
246
247
248#=============== End of Shallow Water Domain ===============================
249
250#Rotation of momentum vector
251def rotate(q, normal, direction = 1):
252    """Rotate the momentum component q (q[1], q[2])
253    from x,y coordinates to coordinates based on normal vector.
254
255    If direction is negative the rotation is inverted.
256
257    Input vector is preserved
258
259    This function is specific to the shallow water wave equation
260    """
261
262    from Numeric import zeros, Float
263
264    assert len(q) == 3,\
265           'Vector of conserved quantities must have length 3'\
266           'for 2D shallow water equation'
267
268    try:
269        l = len(normal)
270    except:
271        raise 'Normal vector must be an Numeric array'
272
273    assert l == 2, 'Normal vector must have 2 components'
274
275
276    n1 = normal[0]
277    n2 = normal[1]
278
279    r = zeros(len(q), Float) #Rotated quantities
280    r[0] = q[0]              #First quantity, height, is not rotated
281
282    if direction == -1:
283        n2 = -n2
284
285
286    r[1] =  n1*q[1] + n2*q[2]
287    r[2] = -n2*q[1] + n1*q[2]
288
289    return r
290
291
292def flux_function(normal, ql, qr, zl, zr):
293    """Compute fluxes between volumes for the shallow water wave equation
294    cast in terms of w = h+z using the 'central scheme' as described in
295
296    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
297    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
298    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
299
300    The implemented formula is given in equation (3.15) on page 714
301
302    Conserved quantities w, uh, are stored as elements 0 and 1
303    in the numerical vectors ql an qr.
304
305    Bed elevations zl and zr.
306    """
307
308    from config import g, epsilon
309    from math import sqrt
310    from Numeric import array
311
312    #print 'ql',ql
313
314    #Align momentums with x-axis
315    #q_left  = rotate(ql, normal, direction = 1)
316    #q_right = rotate(qr, normal, direction = 1)
317    q_left = ql
318    q_left[1] = q_left[1]*normal
319    q_right = qr
320    q_right[1] = q_right[1]*normal
321
322    #z = (zl+zr)/2 #Take average of field values
323    z = 0.5*(zl+zr) #Take average of field values
324
325    w_left  = q_left[0]   #w=h+z
326    h_left  = w_left-z
327    uh_left = q_left[1]
328
329    if h_left < epsilon:
330        u_left = 0.0  #Could have been negative
331        h_left = 0.0
332    else:
333        u_left  = uh_left/h_left
334
335
336    w_right  = q_right[0]  #w=h+z
337    h_right  = w_right-z
338    uh_right = q_right[1]
339
340
341    if h_right < epsilon:
342        u_right = 0.0  #Could have been negative
343        h_right = 0.0
344    else:
345        u_right  = uh_right/h_right
346
347    #vh_left  = q_left[2]
348    #vh_right = q_right[2]
349
350    #print h_right
351    #print u_right
352    #print h_left
353    #print u_right
354
355    soundspeed_left  = sqrt(g*h_left)
356    soundspeed_right = sqrt(g*h_right)
357
358    #Maximal wave speed
359    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
360   
361    #Minimal wave speed
362    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
363   
364    #Flux computation
365
366    #flux_left  = array([u_left*h_left,
367    #                    u_left*uh_left + 0.5*g*h_left**2])
368    #flux_right = array([u_right*h_right,
369    #                    u_right*uh_right + 0.5*g*h_right**2])
370    flux_left  = array([u_left*h_left,
371                        u_left*uh_left + 0.5*g*h_left*h_left])
372    flux_right = array([u_right*h_right,
373                        u_right*uh_right + 0.5*g*h_right*h_right])
374
375    denom = s_max-s_min
376    if denom == 0.0:
377        edgeflux = array([0.0, 0.0])
378        max_speed = 0.0
379    else:
380        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
381        edgeflux += s_max*s_min*(q_right-q_left)/denom
382       
383        edgeflux[1] = edgeflux[1]*normal
384
385        max_speed = max(abs(s_max), abs(s_min))
386
387    return edgeflux, max_speed
388
389
390 
391
392def compute_timestep(domain):
393    import sys
394    from Numeric import zeros, Float
395
396    N = domain.number_of_elements
397
398    #Shortcuts
399    Stage = domain.quantities['stage']
400    Xmom = domain.quantities['xmomentum']
401    Bed = domain.quantities['elevation']
402   
403    stage = Stage.vertex_values
404    xmom =  Xmom.vertex_values
405    bed =   Bed.vertex_values
406
407    stage_bdry = Stage.boundary_values
408    xmom_bdry =  Xmom.boundary_values
409
410    flux = zeros(2, Float) #Work array for summing up fluxes
411    ql = zeros(2, Float)
412    qr = zeros(2, Float)
413
414    #Loop
415    timestep = float(sys.maxint)
416    enter = True
417    for k in range(N):
418
419        flux[:] = 0.  #Reset work array
420        for i in range(2):
421            #Quantities inside volume facing neighbour i
422            ql = [stage[k, i], xmom[k, i]]
423            zl = bed[k, i]
424
425            #Quantities at neighbour on nearest face
426            n = domain.neighbours[k,i]
427            if n < 0:
428                m = -n-1 #Convert negative flag to index
429                qr[0] = stage_bdry[m]
430                qr[1] = xmom_bdry[m]
431                zr = zl #Extend bed elevation to boundary
432            else:
433                #m = domain.neighbour_edges[k,i]
434                m = domain.neighbour_vertices[k,i]
435                qr[0] = stage[n, m]
436                qr[1] =  xmom[n, m]
437                zr = bed[n, m]
438
439
440            #Outward pointing normal vector
441            normal = domain.normals[k, i]
442
443            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
444
445            #Update optimal_timestep
446            try:
447                timestep = min(timestep, domain.cfl*0.5*domain.areas[k]/max_speed)
448            except ZeroDivisionError:
449                pass
450
451    domain.timestep = timestep
452
453def compute_fluxes(domain):
454    """Compute all fluxes and the timestep suitable for all volumes
455    in domain.
456
457    Compute total flux for each conserved quantity using "flux_function"
458
459    Fluxes across each edge are scaled by edgelengths and summed up
460    Resulting flux is then scaled by area and stored in
461    explicit_update for each of the three conserved quantities
462    stage, xmomentum and ymomentum
463
464    The maximal allowable speed computed by the flux_function for each volume
465    is converted to a timestep that must not be exceeded. The minimum of
466    those is computed as the next overall timestep.
467
468    Post conditions:
469      domain.explicit_update is reset to computed flux values
470      domain.timestep is set to the largest step satisfying all volumes.
471    """
472
473    import sys
474    from Numeric import zeros, Float
475
476    N = domain.number_of_elements
477
478    #Shortcuts
479    Stage = domain.quantities['stage']
480    Xmom = domain.quantities['xmomentum']
481#    Ymom = domain.quantities['ymomentum']
482    Bed = domain.quantities['elevation']
483
484    #Arrays
485    #stage = Stage.edge_values
486    #xmom =  Xmom.edge_values
487 #   ymom =  Ymom.edge_values
488    #bed =   Bed.edge_values
489   
490    stage = Stage.vertex_values
491    xmom =  Xmom.vertex_values
492    bed =   Bed.vertex_values
493   
494    #print 'stage edge values', stage
495    #print 'xmom edge values', xmom
496    #print 'bed values', bed
497
498    stage_bdry = Stage.boundary_values
499    xmom_bdry =  Xmom.boundary_values
500    #print 'stage_bdry',stage_bdry
501    #print 'xmom_bdry', xmom_bdry
502#    ymom_bdry =  Ymom.boundary_values
503
504#    flux = zeros(3, Float) #Work array for summing up fluxes
505    flux = zeros(2, Float) #Work array for summing up fluxes
506    ql = zeros(2, Float)
507    qr = zeros(2, Float)
508
509    #Loop
510    timestep = float(sys.maxint)
511    enter = True
512    for k in range(N):
513
514        flux[:] = 0.  #Reset work array
515        #for i in range(3):
516        for i in range(2):
517            #Quantities inside volume facing neighbour i
518            #ql[0] = stage[k, i]
519            #ql[1] = xmom[k, i]
520            ql = [stage[k, i], xmom[k, i]]
521            zl = bed[k, i]
522
523            #Quantities at neighbour on nearest face
524            n = domain.neighbours[k,i]
525            if n < 0:
526                m = -n-1 #Convert negative flag to index
527                qr[0] = stage_bdry[m]
528                qr[1] = xmom_bdry[m]
529                zr = zl #Extend bed elevation to boundary
530            else:
531                #m = domain.neighbour_edges[k,i]
532                m = domain.neighbour_vertices[k,i]
533                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
534                qr[0] = stage[n, m]
535                qr[1] = xmom[n, m]
536                zr = bed[n, m]
537
538
539            #Outward pointing normal vector
540            normal = domain.normals[k, i]
541       
542            #Flux computation using provided function
543            #edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
544            #print 'ql',ql
545            #print 'qr',qr
546           
547
548            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
549
550            #print 'edgeflux', edgeflux
551
552            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
553            # flux = edgefluxleft - edgefluxright
554            flux -= edgeflux #* domain.edgelengths[k,i]
555            #Update optimal_timestep
556            try:
557                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
558                timestep = min(timestep, domain.cfl*0.5*domain.areas[k]/max_speed)
559            except ZeroDivisionError:
560                pass
561
562        #Normalise by area and store for when all conserved
563        #quantities get updated
564        flux /= domain.areas[k]
565
566        Stage.explicit_update[k] = flux[0]
567        Xmom.explicit_update[k] =  flux[1]
568        #Ymom.explicit_update[k] = flux[2]
569        #print "flux cell",k,flux[0]
570
571    domain.timestep = timestep
572    #print domain.quantities['stage'].centroid_values
573
574####################################
575
576# Module functions for gradient limiting (distribute_to_vertices_and_edges)
577
578def distribute_to_vertices_and_edges(domain):
579    """Distribution from centroids to vertices specific to the
580    shallow water wave
581    equation.
582
583    It will ensure that h (w-z) is always non-negative even in the
584    presence of steep bed-slopes by taking a weighted average between shallow
585    and deep cases.
586
587    In addition, all conserved quantities get distributed as per either a
588    constant (order==1) or a piecewise linear function (order==2).
589
590    FIXME: more explanation about removal of artificial variability etc
591
592    Precondition:
593      All quantities defined at centroids and bed elevation defined at
594      vertices.
595
596    Postcondition
597      Conserved quantities defined at vertices
598
599    """
600
601    #from config import optimised_gradient_limiter
602
603    #Remove very thin layers of water
604    protect_against_infinitesimal_and_negative_heights(domain)
605   
606
607    #Extrapolate all conserved quantities
608    #if optimised_gradient_limiter:
609    #    #MH090605 if second order,
610    #    #perform the extrapolation and limiting on
611    #    #all of the conserved quantities
612
613    #    if (domain.order == 1):
614    #        for name in domain.conserved_quantities:
615    #            Q = domain.quantities[name]
616    #            Q.extrapolate_first_order()
617    #    elif domain.order == 2:
618    #        domain.extrapolate_second_order_sw()
619    #    else:
620    #        raise 'Unknown order'
621    #else:
622        #old code:
623
624    for name in domain.conserved_quantities:
625        Q = domain.quantities[name]
626        if domain.order == 1:
627            Q.extrapolate_first_order()
628        elif domain.order == 2:
629            #print "add extrapolate second order to shallow water"
630            #if name != 'height':
631            Q.extrapolate_second_order()
632            #Q.limit()
633        else:
634            raise 'Unknown order'
635
636    #Take bed elevation into account when water heights are small
637    #balance_deep_and_shallow(domain)
638    #protect_against_infinitesimal_and_negative_heights(domain)
639
640#Compute edge values by interpolation
641    #for name in domain.conserved_quantities:
642    #    Q = domain.quantities[name]
643    #    Q.interpolate_from_vertices_to_edges()   
644   
645def protect_against_infinitesimal_and_negative_heights(domain):
646    """Protect against infinitesimal heights and associated high velocities
647    """
648
649    #Shortcuts
650    wc = domain.quantities['stage'].centroid_values
651    zc = domain.quantities['elevation'].centroid_values
652    xmomc = domain.quantities['xmomentum'].centroid_values
653#    ymomc = domain.quantities['ymomentum'].centroid_values
654    hc = wc - zc  #Water depths at centroids
655
656    zv = domain.quantities['elevation'].vertex_values
657    wv = domain.quantities['stage'].vertex_values
658    hv = wv-zv
659    xmomv = domain.quantities['xmomentum'].vertex_values
660    #remove the above two lines and corresponding code below
661
662    #Update
663    for k in range(domain.number_of_elements):
664
665        if hc[k] < domain.minimum_allowed_height:
666            #Control stage
667            if hc[k] < domain.epsilon:
668                wc[k] = zc[k] # Contain 'lost mass' error
669                wv[k,0] = zv[k,0]
670                wv[k,1] = zv[k,1]
671               
672            xmomc[k] = 0.0
673           
674        #N = domain.number_of_elements
675        #if (k == 0) | (k==N-1):
676        #    wc[k] = zc[k] # Contain 'lost mass' error
677        #    wv[k,0] = zv[k,0]
678        #    wv[k,1] = zv[k,1]
679   
680def h_limiter(domain):
681    """Limit slopes for each volume to eliminate artificial variance
682    introduced by e.g. second order extrapolator
683
684    limit on h = w-z
685
686    This limiter depends on two quantities (w,z) so it resides within
687    this module rather than within quantity.py
688    """
689
690    from Numeric import zeros, Float
691
692    N = domain.number_of_elements
693    beta_h = domain.beta_h
694
695    #Shortcuts
696    wc = domain.quantities['stage'].centroid_values
697    zc = domain.quantities['elevation'].centroid_values
698    hc = wc - zc
699
700    wv = domain.quantities['stage'].vertex_values
701    zv = domain.quantities['elevation'].vertex_values
702    hv = wv-zv
703
704    hvbar = zeros(hv.shape, Float) #h-limited values
705
706    #Find min and max of this and neighbour's centroid values
707    hmax = zeros(hc.shape, Float)
708    hmin = zeros(hc.shape, Float)
709
710    for k in range(N):
711        hmax[k] = hmin[k] = hc[k]
712        #for i in range(3):
713        for i in range(2):   
714            n = domain.neighbours[k,i]
715            if n >= 0:
716                hn = hc[n] #Neighbour's centroid value
717
718                hmin[k] = min(hmin[k], hn)
719                hmax[k] = max(hmax[k], hn)
720
721
722    #Diffences between centroids and maxima/minima
723    dhmax = hmax - hc
724    dhmin = hmin - hc
725
726    #Deltas between vertex and centroid values
727    dh = zeros(hv.shape, Float)
728    #for i in range(3):
729    for i in range(2):
730        dh[:,i] = hv[:,i] - hc
731
732    #Phi limiter
733    for k in range(N):
734
735        #Find the gradient limiter (phi) across vertices
736        phi = 1.0
737        #for i in range(3):
738        for i in range(2):
739            r = 1.0
740            if (dh[k,i] > 0): r = dhmax[k]/dh[k,i]
741            if (dh[k,i] < 0): r = dhmin[k]/dh[k,i]
742
743            phi = min( min(r*beta_h, 1), phi )
744
745        #Then update using phi limiter
746        #for i in range(3):
747        for i in range(2):
748            hvbar[k,i] = hc[k] + phi*dh[k,i]
749
750    return hvbar
751
752def balance_deep_and_shallow(domain):
753    """Compute linear combination between stage as computed by
754    gradient-limiters limiting using w, and stage computed by
755    gradient-limiters limiting using h (h-limiter).
756    The former takes precedence when heights are large compared to the
757    bed slope while the latter takes precedence when heights are
758    relatively small.  Anything in between is computed as a balanced
759    linear combination in order to avoid numerical disturbances which
760    would otherwise appear as a result of hard switching between
761    modes.
762
763    The h-limiter is always applied irrespective of the order.
764    """
765
766    #Shortcuts
767    wc = domain.quantities['stage'].centroid_values
768    zc = domain.quantities['elevation'].centroid_values
769    hc = wc - zc
770
771    wv = domain.quantities['stage'].vertex_values
772    zv = domain.quantities['elevation'].vertex_values
773    hv = wv-zv
774
775    #Limit h
776    hvbar = h_limiter(domain)
777
778    for k in range(domain.number_of_elements):
779        #Compute maximal variation in bed elevation
780        #  This quantitiy is
781        #    dz = max_i abs(z_i - z_c)
782        #  and it is independent of dimension
783        #  In the 1d case zc = (z0+z1)/2
784        #  In the 2d case zc = (z0+z1+z2)/3
785
786        dz = max(abs(zv[k,0]-zc[k]),
787                 abs(zv[k,1]-zc[k]))#,
788#                 abs(zv[k,2]-zc[k]))
789
790
791        hmin = min( hv[k,:] )
792
793        #Create alpha in [0,1], where alpha==0 means using the h-limited
794        #stage and alpha==1 means using the w-limited stage as
795        #computed by the gradient limiter (both 1st or 2nd order)
796
797        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
798        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
799
800        if dz > 0.0:
801            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
802        else:
803            #Flat bed
804            alpha = 1.0
805
806        alpha  = 0.0
807        #Let
808        #
809        #  wvi be the w-limited stage (wvi = zvi + hvi)
810        #  wvi- be the h-limited state (wvi- = zvi + hvi-)
811        #
812        #
813        #where i=0,1,2 denotes the vertex ids
814        #
815        #Weighted balance between w-limited and h-limited stage is
816        #
817        #  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
818        #
819        #It follows that the updated wvi is
820        #  wvi := zvi + (1-alpha)*hvi- + alpha*hvi
821        #
822        # Momentum is balanced between constant and limited
823
824
825        #for i in range(3):
826        #    wv[k,i] = zv[k,i] + hvbar[k,i]
827
828        #return
829
830        if alpha < 1:
831
832            #for i in range(3):
833            for i in range(2):
834                wv[k,i] = zv[k,i] + (1.0-alpha)*hvbar[k,i] + alpha*hv[k,i]
835
836            #Momentums at centroids
837            xmomc = domain.quantities['xmomentum'].centroid_values
838#            ymomc = domain.quantities['ymomentum'].centroid_values
839
840            #Momentums at vertices
841            xmomv = domain.quantities['xmomentum'].vertex_values
842#           ymomv = domain.quantities['ymomentum'].vertex_values
843
844            # Update momentum as a linear combination of
845            # xmomc and ymomc (shallow) and momentum
846            # from extrapolator xmomv and ymomv (deep).
847            xmomv[k,:] = (1.0-alpha)*xmomc[k] + alpha*xmomv[k,:]
848#            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
849
850
851###############################################
852#Boundaries - specific to the shallow water wave equation
853class Reflective_boundary(Boundary):
854    """Reflective boundary returns same conserved quantities as
855    those present in its neighbour volume but reflected.
856
857    This class is specific to the shallow water equation as it
858    works with the momentum quantities assumed to be the second
859    and third conserved quantities.
860    """
861
862    def __init__(self, domain = None):
863        Boundary.__init__(self)
864
865        if domain is None:
866            msg = 'Domain must be specified for reflective boundary'
867            raise msg
868
869        #Handy shorthands
870        #self.stage   = domain.quantities['stage'].edge_values
871        #self.xmom    = domain.quantities['xmomentum'].edge_values
872        #self.ymom    = domain.quantities['ymomentum'].edge_values
873        self.normals = domain.normals
874        self.stage   = domain.quantities['stage'].vertex_values
875        self.xmom    = domain.quantities['xmomentum'].vertex_values
876
877        from Numeric import zeros, Float
878        #self.conserved_quantities = zeros(3, Float)
879        self.conserved_quantities = zeros(2, Float)
880
881    def __repr__(self):
882        return 'Reflective_boundary'
883
884
885    def evaluate(self, vol_id, edge_id):
886        """Reflective boundaries reverses the outward momentum
887        of the volume they serve.
888        """
889
890        q = self.conserved_quantities
891        q[0] = self.stage[vol_id, edge_id]
892        q[1] = self.xmom[vol_id, edge_id]
893        #q[2] = self.ymom[vol_id, edge_id]
894        #normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
895        #normal = self.normals[vol_id, 2*edge_id:2*edge_id+1]
896        normal = self.normals[vol_id,edge_id]
897
898        #r = rotate(q, normal, direction = 1)
899        #r[1] = -r[1]
900        #q = rotate(r, normal, direction = -1)
901        r = q
902        r[1] = normal*r[1]
903        r[1] = -r[1]
904        r[1] = normal*r[1]
905        q = r
906        #For start interval there is no outward momentum so do not need to
907        #reverse direction in this case
908
909        return q
910
911class Dirichlet_boundary(Boundary):
912    """Dirichlet boundary returns constant values for the
913    conserved quantities
914    """
915
916
917    def __init__(self, conserved_quantities=None):
918        Boundary.__init__(self)
919
920        if conserved_quantities is None:
921            msg = 'Must specify one value for each conserved quantity'
922            raise msg
923
924        from Numeric import array, Float
925        self.conserved_quantities=array(conserved_quantities).astype(Float)
926
927    def __repr__(self):
928        return 'Dirichlet boundary (%s)' %self.conserved_quantities
929
930    def evaluate(self, vol_id=None, edge_id=None):
931        return self.conserved_quantities
932
933
934#########################
935#Standard forcing terms:
936#
937def gravity(domain):
938    """Apply gravitational pull in the presence of bed slope
939    """
940
941    from util import gradient
942    from Numeric import zeros, Float, array, sum
943
944    xmom = domain.quantities['xmomentum'].explicit_update
945    stage = domain.quantities['stage'].explicit_update
946#    ymom = domain.quantities['ymomentum'].explicit_update
947
948    Stage = domain.quantities['stage']
949    Elevation = domain.quantities['elevation']
950    #h = Stage.edge_values - Elevation.edge_values
951    h = Stage.vertex_values - Elevation.vertex_values
952    b = Elevation.vertex_values
953    w = Stage.vertex_values
954
955    x = domain.get_vertex_coordinates()
956    g = domain.g
957
958    for k in range(domain.number_of_elements):
959#        avg_h = sum( h[k,:] )/3
960        avg_h = sum( h[k,:] )/2
961
962        #Compute bed slope
963        #x0, y0, x1, y1, x2, y2 = x[k,:]
964        x0, x1 = x[k,:]
965        #z0, z1, z2 = v[k,:]
966        b0, b1 = b[k,:]
967
968        w0, w1 = w[k,:]
969        wx = gradient(x0, x1, w0, w1)
970
971        #zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
972        bx = gradient(x0, x1, b0, b1)
973       
974        #Update momentum (explicit update is reset to source values)
975        xmom[k] += -g*bx*avg_h
976        #xmom[k] = -g*bx*avg_h
977        #stage[k] = 0.0
978 
979 
980def manning_friction(domain):
981    """Apply (Manning) friction to water momentum
982    """
983
984    from math import sqrt
985
986    w = domain.quantities['stage'].centroid_values
987    z = domain.quantities['elevation'].centroid_values
988    h = w-z
989
990    uh = domain.quantities['xmomentum'].centroid_values
991    #vh = domain.quantities['ymomentum'].centroid_values
992    eta = domain.quantities['friction'].centroid_values
993
994    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
995    #ymom_update = domain.quantities['ymomentum'].semi_implicit_update
996
997    N = domain.number_of_elements
998    eps = domain.minimum_allowed_height
999    g = domain.g
1000
1001    for k in range(N):
1002        if eta[k] >= eps:
1003            if h[k] >= eps:
1004                #S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
1005                S = -g * eta[k]**2 * uh[k]
1006                S /= h[k]**(7.0/3)
1007
1008                #Update momentum
1009                xmom_update[k] += S*uh[k]
1010                #ymom_update[k] += S*vh[k]
1011
1012def linear_friction(domain):
1013    """Apply linear friction to water momentum
1014
1015    Assumes quantity: 'linear_friction' to be present
1016    """
1017
1018    from math import sqrt
1019
1020    w = domain.quantities['stage'].centroid_values
1021    z = domain.quantities['elevation'].centroid_values
1022    h = w-z
1023
1024    uh = domain.quantities['xmomentum'].centroid_values
1025#    vh = domain.quantities['ymomentum'].centroid_values
1026    tau = domain.quantities['linear_friction'].centroid_values
1027
1028    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1029#    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1030
1031    N = domain.number_of_elements
1032    eps = domain.minimum_allowed_height
1033    g = domain.g #Not necessary? Why was this added?
1034
1035    for k in range(N):
1036        if tau[k] >= eps:
1037            if h[k] >= eps:
1038                S = -tau[k]/h[k]
1039
1040                #Update momentum
1041                xmom_update[k] += S*uh[k]
1042 #              ymom_update[k] += S*vh[k]
1043
1044
1045
1046def check_forcefield(f):
1047    """Check that f is either
1048    1: a callable object f(t,x,y), where x and y are vectors
1049       and that it returns an array or a list of same length
1050       as x and y
1051    2: a scalar
1052    """
1053
1054    from Numeric import ones, Float, array
1055
1056
1057    if callable(f):
1058        #N = 3
1059        N = 2
1060        #x = ones(3, Float)
1061        #y = ones(3, Float)
1062        x = ones(2, Float)
1063        #y = ones(2, Float)
1064       
1065        try:
1066            #q = f(1.0, x=x, y=y)
1067            q = f(1.0, x=x)
1068        except Exception, e:
1069            msg = 'Function %s could not be executed:\n%s' %(f, e)
1070            #FIXME: Reconsider this semantics
1071            raise msg
1072
1073        try:
1074            q = array(q).astype(Float)
1075        except:
1076            msg = 'Return value from vector function %s could ' %f
1077            msg += 'not be converted into a Numeric array of floats.\n'
1078            msg += 'Specified function should return either list or array.'
1079            raise msg
1080
1081        #Is this really what we want?
1082        msg = 'Return vector from function %s ' %f
1083        msg += 'must have same lenght as input vectors'
1084        assert len(q) == N, msg
1085
1086    else:
1087        try:
1088            f = float(f)
1089        except:
1090            msg = 'Force field %s must be either a scalar' %f
1091            msg += ' or a vector function'
1092            raise msg
1093    return f
1094
1095class Wind_stress:
1096    """Apply wind stress to water momentum in terms of
1097    wind speed [m/s] and wind direction [degrees]
1098    """
1099
1100    def __init__(self, *args, **kwargs):
1101        """Initialise windfield from wind speed s [m/s]
1102        and wind direction phi [degrees]
1103
1104        Inputs v and phi can be either scalars or Python functions, e.g.
1105
1106        W = Wind_stress(10, 178)
1107
1108        #FIXME - 'normal' degrees are assumed for now, i.e. the
1109        vector (1,0) has zero degrees.
1110        We may need to convert from 'compass' degrees later on and also
1111        map from True north to grid north.
1112
1113        Arguments can also be Python functions of t,x,y as in
1114
1115        def speed(t,x,y):
1116            ...
1117            return s
1118
1119        def angle(t,x,y):
1120            ...
1121            return phi
1122
1123        where x and y are vectors.
1124
1125        and then pass the functions in
1126
1127        W = Wind_stress(speed, angle)
1128
1129        The instantiated object W can be appended to the list of
1130        forcing_terms as in
1131
1132        Alternatively, one vector valued function for (speed, angle)
1133        can be applied, providing both quantities simultaneously.
1134        As in
1135        W = Wind_stress(F), where returns (speed, angle) for each t.
1136
1137        domain.forcing_terms.append(W)
1138        """
1139
1140        from config import rho_a, rho_w, eta_w
1141        from Numeric import array, Float
1142
1143        if len(args) == 2:
1144            s = args[0]
1145            phi = args[1]
1146        elif len(args) == 1:
1147            #Assume vector function returning (s, phi)(t,x,y)
1148            vector_function = args[0]
1149            #s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
1150            #phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
1151            s = lambda t,x: vector_function(t,x=x)[0]
1152            phi = lambda t,x: vector_function(t,x=x)[1]
1153        else:
1154           #Assume info is in 2 keyword arguments
1155
1156           if len(kwargs) == 2:
1157               s = kwargs['s']
1158               phi = kwargs['phi']
1159           else:
1160               raise 'Assumes two keyword arguments: s=..., phi=....'
1161
1162        print 'phi', phi
1163        self.speed = check_forcefield(s)
1164        self.phi = check_forcefield(phi)
1165
1166        self.const = eta_w*rho_a/rho_w
1167
1168
1169    def __call__(self, domain):
1170        """Evaluate windfield based on values found in domain
1171        """
1172
1173        from math import pi, cos, sin, sqrt
1174        from Numeric import Float, ones, ArrayType
1175
1176        xmom_update = domain.quantities['xmomentum'].explicit_update
1177        #ymom_update = domain.quantities['ymomentum'].explicit_update
1178
1179        N = domain.number_of_elements
1180        t = domain.time
1181
1182        if callable(self.speed):
1183            xc = domain.get_centroid_coordinates()
1184            #s_vec = self.speed(t, xc[:,0], xc[:,1])
1185            s_vec = self.speed(t, xc)
1186        else:
1187            #Assume s is a scalar
1188
1189            try:
1190                s_vec = self.speed * ones(N, Float)
1191            except:
1192                msg = 'Speed must be either callable or a scalar: %s' %self.s
1193                raise msg
1194
1195
1196        if callable(self.phi):
1197            xc = domain.get_centroid_coordinates()
1198            #phi_vec = self.phi(t, xc[:,0], xc[:,1])
1199            phi_vec = self.phi(t, xc)
1200        else:
1201            #Assume phi is a scalar
1202
1203            try:
1204                phi_vec = self.phi * ones(N, Float)
1205            except:
1206                msg = 'Angle must be either callable or a scalar: %s' %self.phi
1207                raise msg
1208
1209        #assign_windfield_values(xmom_update, ymom_update,
1210        #                        s_vec, phi_vec, self.const)
1211        assign_windfield_values(xmom_update, s_vec, phi_vec, self.const)
1212
1213
1214#def assign_windfield_values(xmom_update, ymom_update,
1215#                            s_vec, phi_vec, const):
1216def assign_windfield_values(xmom_update, s_vec, phi_vec, const):
1217    """Python version of assigning wind field to update vectors.
1218    A c version also exists (for speed)
1219    """
1220    from math import pi, cos, sin, sqrt
1221
1222    N = len(s_vec)
1223    for k in range(N):
1224        s = s_vec[k]
1225        phi = phi_vec[k]
1226
1227        #Convert to radians
1228        phi = phi*pi/180
1229
1230        #Compute velocity vector (u, v)
1231        u = s*cos(phi)
1232        v = s*sin(phi)
1233
1234        #Compute wind stress
1235        #S = const * sqrt(u**2 + v**2)
1236        S = const * u
1237        xmom_update[k] += S*u
1238        #ymom_update[k] += S*v
Note: See TracBrowser for help on using the repository browser.