source: inundation/pyvolution/shallow_water.py @ 2566

Last change on this file since 2566 was 2566, checked in by ole, 18 years ago

Created maximum_allowed_speed in config and passed it through.
Also ditched epsilon and used h<=0 instead.

  • Property svn:keywords set to LastChangedDate LastChangedRevision LastChangedBy HeadURL Id
File size: 53.8 KB
RevLine 
[2566]1"""Class Domain -
22D triangular 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 + G_y = S
10
11where
12
13U = [w, uh, vh]
14E = [uh, u^2h + gh^2/2, uvh]
15G = [vh, uvh, v^2h + gh^2/2]
16S represents source terms forcing the system
17(e.g. gravity, friction, wind stress, ...)
18
19and _t, _x, _y denote the derivative with respect to t, x and y respectively.
20
21The quantities are
22
23symbol    variable name    explanation
24x         x                horizontal distance from origin [m]
25y         y                vertical distance from origin [m]
26z         elevation        elevation of bed on which flow is modelled [m]
27h         height           water height above z [m]
28w         stage            absolute water level, w = z+h [m]
29u                          speed in the x direction [m/s]
30v                          speed in the y direction [m/s]
31uh        xmomentum        momentum in the x direction [m^2/s]
32vh        ymomentum        momentum in the y direction [m^2/s]
33
34eta                        mannings friction coefficient [to appear]
35nu                         wind stress coefficient [to appear]
36
37The conserved quantities are w, uh, vh
38
39
40For details see e.g.
41Christopher Zoppou and Stephen Roberts,
42Catastrophic Collapse of Water Supply Reservoirs in Urban Areas,
43Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999
44
45
46Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
47Geoscience Australia, 2004
48"""
49
50#Subversion keywords:
51#
52#$LastChangedDate: 2006-03-21 03:25:40 +0000 (Tue, 21 Mar 2006) $
53#$LastChangedRevision: 2566 $
54#$LastChangedBy: ole $
55
56
57from domain import Domain as Generic_Domain
58from generic_boundary_conditions import Boundary
59from generic_boundary_conditions import File_boundary
60from generic_boundary_conditions import Dirichlet_boundary
61from generic_boundary_conditions import Time_boundary
62from generic_boundary_conditions import Transmissive_boundary
63
64from utilities.numerical_tools import gradient, mean
65
66
67
68#Shallow water domain
69class Domain(Generic_Domain):
70
71    def __init__(self, coordinates, vertices, boundary = None,
72                 tagged_elements = None, geo_reference = None,
73                 use_inscribed_circle=False):
74
75        conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
76        other_quantities = ['elevation', 'friction']
77        Generic_Domain.__init__(self, coordinates, vertices, boundary,
78                                conserved_quantities, other_quantities,
79                                tagged_elements, geo_reference, use_inscribed_circle)
80
81        from config import minimum_allowed_height, maximum_allowed_speed, g
82        self.minimum_allowed_height = minimum_allowed_height
83        self.maximum_allowed_speed = maximum_allowed_speed     
84        self.g = g
85
86        self.forcing_terms.append(gravity)
87        self.forcing_terms.append(manning_friction)
88
89        #Realtime visualisation
90        self.visualiser = None
91        self.visualise  = False
92        self.visualise_color_stage = False
93        self.visualise_stage_range = 1.0
94        self.visualise_timer = True
95        self.visualise_range_z = None
96
97        #Stored output
98        self.store = True
99        self.format = 'sww'
100        self.smooth = True
101
102        #Reduction operation for get_vertex_values
103        self.reduction = mean
104        #self.reduction = min  #Looks better near steep slopes
105
106        self.quantities_to_be_stored = ['stage','xmomentum','ymomentum']
107
108
109    def set_quantities_to_be_stored(self, q):
110        """Specify which quantities will be stored in the sww file.
111
112        q must be either:
113          - the name of a quantity
114          - a list of quantity names
115          - None
116
117        In the two first cases, the named quantities will be stored at each yieldstep
118        (This is in addition to the quantities elevation and friction) 
119        If q is None, storage will be switched off altogether.
120        """
121
122
123        if q is None:
124            self.quantities_to_be_stored = []           
125            self.store = False
126            return
127
128        if isinstance(q, basestring):
129            q = [q] # Turn argument into a list
130
131        #Check correcness   
132        for quantity_name in q:
133            msg = 'Quantity %s is not a valid conserved quantity' %quantity_name
134            assert quantity_name in self.conserved_quantities, msg
135       
136        self.quantities_to_be_stored = q
137       
138
139    def initialise_visualiser(self,scale_z=1.0,rect=None):
140        #Realtime visualisation
141        if self.visualiser is None:
142            from realtime_visualisation_new import Visualiser
143            self.visualiser = Visualiser(self,scale_z,rect)
144            self.visualiser.setup['elevation']=True
145            self.visualiser.updating['stage']=True
146        self.visualise = True
147        if self.visualise_color_stage == True:
148            self.visualiser.coloring['stage'] = True
149            self.visualiser.qcolor['stage'] = (0.0, 0.0, 0.8)
150        print 'initialise visualiser'
151        print self.visualiser.setup
152        print self.visualiser.updating
153
154    def check_integrity(self):
155        Generic_Domain.check_integrity(self)
156
157        #Check that we are solving the shallow water wave equation
158
159        msg = 'First conserved quantity must be "stage"'
160        assert self.conserved_quantities[0] == 'stage', msg
161        msg = 'Second conserved quantity must be "xmomentum"'
162        assert self.conserved_quantities[1] == 'xmomentum', msg
163        msg = 'Third conserved quantity must be "ymomentum"'
164        assert self.conserved_quantities[2] == 'ymomentum', msg
165
166    def extrapolate_second_order_sw(self):
167        #Call correct module function
168        #(either from this module or C-extension)
169        extrapolate_second_order_sw(self)
170
171    def compute_fluxes(self):
172        #Call correct module function
173        #(either from this module or C-extension)
174        compute_fluxes(self)
175
176    def distribute_to_vertices_and_edges(self):
177        #Call correct module function
178        #(either from this module or C-extension)
179        distribute_to_vertices_and_edges(self)
180
181
182    #FIXME: Under construction
183#     def set_defaults(self):
184#         """Set default values for uninitialised quantities.
185#         This is specific to the shallow water wave equation
186#         Defaults for 'elevation', 'friction', 'xmomentum' and 'ymomentum'
187#         are 0.0. Default for 'stage' is whatever the value of 'elevation'.
188#         """
189
190#         for name in self.other_quantities + self.conserved_quantities:
191#             print name
192#             print self.quantities.keys()
193#             if not self.quantities.has_key(name):
194#                 if name == 'stage':
195
196#                     if self.quantities.has_key('elevation'):
197#                         z = self.quantities['elevation'].vertex_values
198#                         self.set_quantity(name, z)
199#                     else:
200#                         self.set_quantity(name, 0.0)
201#                 else:
202#                     self.set_quantity(name, 0.0)
203
204
205
206#         #Lift negative heights up
207#         #z = self.quantities['elevation'].vertex_values
208#         #w = self.quantities['stage'].vertex_values
209
210#         #h = w-z
211
212#         #for k in range(h.shape[0]):
213#         #    for i in range(3):
214#         #        if h[k, i] < 0.0:
215#         #            w[k, i] = z[k, i]
216
217
218#         #self.quantities['stage'].interpolate()
219
220
221    def evolve(self,
222               yieldstep = None,
223               finaltime = None,
224               duration = None,               
225               skip_initial_step = False):
226        """Specialisation of basic evolve method from parent class
227        """
228
229        #Call check integrity here rather than from user scripts
230        #self.check_integrity()
231
232        msg = 'Parameter beta_h must be in the interval [0, 1['
233        assert 0 <= self.beta_h < 1.0, msg
234        msg = 'Parameter beta_w must be in the interval [0, 1['
235        assert 0 <= self.beta_w < 1.0, msg
236
237
238        #Initial update of vertex and edge values before any storage
239        #and or visualisation
240        self.distribute_to_vertices_and_edges()
241
242        #Initialise real time viz if requested
243        if self.visualise is True and self.time == 0.0:
244            if self.visualiser is None:
245                self.initialise_visualiser()
246
247            self.visualiser.update_timer()
248            self.visualiser.setup_all()
249
250        #Store model data, e.g. for visualisation
251        if self.store is True and self.time == 0.0:
252            self.initialise_storage()
253            #print 'Storing results in ' + self.writer.filename
254        else:
255            pass
256            #print 'Results will not be stored.'
257            #print 'To store results set domain.store = True'
258            #FIXME: Diagnostic output should be controlled by
259            # a 'verbose' flag living in domain (or in a parent class)
260
261        #Call basic machinery from parent class
262        for t in Generic_Domain.evolve(self,
263                                       yieldstep=yieldstep,
264                                       finaltime=finaltime,
265                                       duration=duration,
266                                       skip_initial_step=skip_initial_step):
267            #Real time viz
268            if self.visualise is True:
269                self.visualiser.update_all()
270                self.visualiser.update_timer()
271
272
273            #Store model data, e.g. for subsequent visualisation
274            if self.store is True:
275                self.store_timestep(self.quantities_to_be_stored)
276
277            #FIXME: Could maybe be taken from specified list
278            #of 'store every step' quantities
279
280            #Pass control on to outer loop for more specific actions
281            yield(t)
282
283    def initialise_storage(self):
284        """Create and initialise self.writer object for storing data.
285        Also, save x,y and bed elevation
286        """
287
288        import data_manager
289
290        #Initialise writer
291        self.writer = data_manager.get_dataobject(self, mode = 'w')
292
293        #Store vertices and connectivity
294        self.writer.store_connectivity()
295
296
297    def store_timestep(self, name):
298        """Store named quantity and time.
299
300        Precondition:
301           self.write has been initialised
302        """
303        self.writer.store_timestep(name)
304
305
306#=============== End of Shallow Water Domain ===============================
307
308
309
310#Rotation of momentum vector
311def rotate(q, normal, direction = 1):
312    """Rotate the momentum component q (q[1], q[2])
313    from x,y coordinates to coordinates based on normal vector.
314
315    If direction is negative the rotation is inverted.
316
317    Input vector is preserved
318
319    This function is specific to the shallow water wave equation
320    """
321
322    from Numeric import zeros, Float
323
324    assert len(q) == 3,\
325           'Vector of conserved quantities must have length 3'\
326           'for 2D shallow water equation'
327
328    try:
329        l = len(normal)
330    except:
331        raise 'Normal vector must be an Numeric array'
332
333    assert l == 2, 'Normal vector must have 2 components'
334
335
336    n1 = normal[0]
337    n2 = normal[1]
338
339    r = zeros(len(q), Float) #Rotated quantities
340    r[0] = q[0]              #First quantity, height, is not rotated
341
342    if direction == -1:
343        n2 = -n2
344
345
346    r[1] =  n1*q[1] + n2*q[2]
347    r[2] = -n2*q[1] + n1*q[2]
348
349    return r
350
351
352
353####################################
354# Flux computation
355def flux_function(normal, ql, qr, zl, zr):
356    """Compute fluxes between volumes for the shallow water wave equation
357    cast in terms of w = h+z using the 'central scheme' as described in
358
359    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
360    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
361    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
362
363    The implemented formula is given in equation (3.15) on page 714
364
365    Conserved quantities w, uh, vh are stored as elements 0, 1 and 2
366    in the numerical vectors ql an qr.
367
368    Bed elevations zl and zr.
369    """
370
371    from config import g, epsilon
372    from math import sqrt
373    from Numeric import array
374
375    #Align momentums with x-axis
376    q_left  = rotate(ql, normal, direction = 1)
377    q_right = rotate(qr, normal, direction = 1)
378
379    z = (zl+zr)/2 #Take average of field values
380
381    w_left  = q_left[0]   #w=h+z
382    h_left  = w_left-z
383    uh_left = q_left[1]
384
385    if h_left < epsilon:
386        u_left = 0.0  #Could have been negative
387        h_left = 0.0
388    else:
389        u_left  = uh_left/h_left
390
391
392    w_right  = q_right[0]  #w=h+z
393    h_right  = w_right-z
394    uh_right = q_right[1]
395
396
397    if h_right < epsilon:
398        u_right = 0.0  #Could have been negative
399        h_right = 0.0
400    else:
401        u_right  = uh_right/h_right
402
403    vh_left  = q_left[2]
404    vh_right = q_right[2]
405
406    soundspeed_left  = sqrt(g*h_left)
407    soundspeed_right = sqrt(g*h_right)
408
409    #Maximal wave speed
410    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
411
412    #Minimal wave speed
413    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
414
415    #Flux computation
416
417    #FIXME(Ole): Why is it again that we don't
418    #use uh_left and uh_right directly in the first entries?
419    flux_left  = array([u_left*h_left,
420                        u_left*uh_left + 0.5*g*h_left**2,
421                        u_left*vh_left])
422    flux_right = array([u_right*h_right,
423                        u_right*uh_right + 0.5*g*h_right**2,
424                        u_right*vh_right])
425
426    denom = s_max-s_min
427    if denom == 0.0:
428        edgeflux = array([0.0, 0.0, 0.0])
429        max_speed = 0.0
430    else:
431        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
432        edgeflux += s_max*s_min*(q_right-q_left)/denom
433
434        edgeflux = rotate(edgeflux, normal, direction=-1)
435        max_speed = max(abs(s_max), abs(s_min))
436
437    return edgeflux, max_speed
438
439
440def compute_fluxes(domain):
441    """Compute all fluxes and the timestep suitable for all volumes
442    in domain.
443
444    Compute total flux for each conserved quantity using "flux_function"
445
446    Fluxes across each edge are scaled by edgelengths and summed up
447    Resulting flux is then scaled by area and stored in
448    explicit_update for each of the three conserved quantities
449    stage, xmomentum and ymomentum
450
451    The maximal allowable speed computed by the flux_function for each volume
452    is converted to a timestep that must not be exceeded. The minimum of
453    those is computed as the next overall timestep.
454
455    Post conditions:
456      domain.explicit_update is reset to computed flux values
457      domain.timestep is set to the largest step satisfying all volumes.
458    """
459
460    import sys
461    from Numeric import zeros, Float
462
463    N = domain.number_of_elements
464
465    #Shortcuts
466    Stage = domain.quantities['stage']
467    Xmom = domain.quantities['xmomentum']
468    Ymom = domain.quantities['ymomentum']
469    Bed = domain.quantities['elevation']
470
471    #Arrays
472    stage = Stage.edge_values
473    xmom =  Xmom.edge_values
474    ymom =  Ymom.edge_values
475    bed =   Bed.edge_values
476
477    stage_bdry = Stage.boundary_values
478    xmom_bdry =  Xmom.boundary_values
479    ymom_bdry =  Ymom.boundary_values
480
481    flux = zeros(3, Float) #Work array for summing up fluxes
482
483
484    #Loop
485    timestep = float(sys.maxint)
486    for k in range(N):
487
488        flux[:] = 0.  #Reset work array
489        for i in range(3):
490            #Quantities inside volume facing neighbour i
491            ql = [stage[k, i], xmom[k, i], ymom[k, i]]
492            zl = bed[k, i]
493
494            #Quantities at neighbour on nearest face
495            n = domain.neighbours[k,i]
496            if n < 0:
497                m = -n-1 #Convert negative flag to index
498                qr = [stage_bdry[m], xmom_bdry[m], ymom_bdry[m]]
499                zr = zl #Extend bed elevation to boundary
500            else:
501                m = domain.neighbour_edges[k,i]
502                qr = [stage[n, m], xmom[n, m], ymom[n, m]]
503                zr = bed[n, m]
504
505
506            #Outward pointing normal vector
507            normal = domain.normals[k, 2*i:2*i+2]
508
509            #Flux computation using provided function
510            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
511            flux -= edgeflux * domain.edgelengths[k,i]
512
513            #Update optimal_timestep
514            try:
515                timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
516            except ZeroDivisionError:
517                pass
518
519        #Normalise by area and store for when all conserved
520        #quantities get updated
521        flux /= domain.areas[k]
522        Stage.explicit_update[k] = flux[0]
523        Xmom.explicit_update[k] = flux[1]
524        Ymom.explicit_update[k] = flux[2]
525
526
527    domain.timestep = timestep
528
529#MH090605 The following method belongs to the shallow_water domain class
530#see comments in the corresponding method in shallow_water_ext.c
531def extrapolate_second_order_sw_c(domain):
532    """Wrapper calling C version of extrapolate_second_order_sw
533    """
534    import sys
535    from Numeric import zeros, Float
536
537    N = domain.number_of_elements
538
539    #Shortcuts
540    Stage = domain.quantities['stage']
541    Xmom = domain.quantities['xmomentum']
542    Ymom = domain.quantities['ymomentum']
543    from shallow_water_ext import extrapolate_second_order_sw
544    extrapolate_second_order_sw(domain,domain.surrogate_neighbours,
545                                domain.number_of_boundaries,
546                                domain.centroid_coordinates,
547                                Stage.centroid_values,
548                                Xmom.centroid_values,
549                                Ymom.centroid_values,
550                                domain.vertex_coordinates,
551                                Stage.vertex_values,
552                                Xmom.vertex_values,
553                                Ymom.vertex_values)
554
555def compute_fluxes_c(domain):
556    """Wrapper calling C version of compute fluxes
557    """
558
559    import sys
560    from Numeric import zeros, Float
561
562    N = domain.number_of_elements
563
564    #Shortcuts
565    Stage = domain.quantities['stage']
566    Xmom = domain.quantities['xmomentum']
567    Ymom = domain.quantities['ymomentum']
568    Bed = domain.quantities['elevation']
569
570    timestep = float(sys.maxint)
571    from shallow_water_ext import compute_fluxes
572    domain.timestep = compute_fluxes(timestep, domain.epsilon, domain.g,
573                                     domain.neighbours,
574                                     domain.neighbour_edges,
575                                     domain.normals,
576                                     domain.edgelengths,
577                                     domain.radii,
578                                     domain.areas,
579                                     Stage.edge_values,
580                                     Xmom.edge_values,
581                                     Ymom.edge_values,
582                                     Bed.edge_values,
583                                     Stage.boundary_values,
584                                     Xmom.boundary_values,
585                                     Ymom.boundary_values,
586                                     Stage.explicit_update,
587                                     Xmom.explicit_update,
588                                     Ymom.explicit_update,
589                                     domain.already_computed_flux)
590
591
592####################################
593# Module functions for gradient limiting (distribute_to_vertices_and_edges)
594
595def distribute_to_vertices_and_edges(domain):
596    """Distribution from centroids to vertices specific to the
597    shallow water wave
598    equation.
599
600    It will ensure that h (w-z) is always non-negative even in the
601    presence of steep bed-slopes by taking a weighted average between shallow
602    and deep cases.
603
604    In addition, all conserved quantities get distributed as per either a
605    constant (order==1) or a piecewise linear function (order==2).
606
607    FIXME: more explanation about removal of artificial variability etc
608
609    Precondition:
610      All quantities defined at centroids and bed elevation defined at
611      vertices.
612
613    Postcondition
614      Conserved quantities defined at vertices
615
616    """
617
618    from config import optimised_gradient_limiter
619
620    #Remove very thin layers of water
621    protect_against_infinitesimal_and_negative_heights(domain)
622
623    #Extrapolate all conserved quantities
624    if optimised_gradient_limiter:
625        #MH090605 if second order,
626        #perform the extrapolation and limiting on
627        #all of the conserved quantitie
628
629        if (domain.order == 1):
630            for name in domain.conserved_quantities:
631                Q = domain.quantities[name]
632                Q.extrapolate_first_order()
633        elif domain.order == 2:
634            domain.extrapolate_second_order_sw()
635        else:
636            raise 'Unknown order'
637    else:
638        #old code:
639        for name in domain.conserved_quantities:
640            Q = domain.quantities[name]
641            if domain.order == 1:
642                Q.extrapolate_first_order()
643            elif domain.order == 2:
644                Q.extrapolate_second_order()
645                Q.limit()
646            else:
647                raise 'Unknown order'
648
649
650    #Take bed elevation into account when water heights are small
651    balance_deep_and_shallow(domain)
652
653    #Compute edge values by interpolation
654    for name in domain.conserved_quantities:
655        Q = domain.quantities[name]
656        Q.interpolate_from_vertices_to_edges()
657
658
659def protect_against_infinitesimal_and_negative_heights(domain):
660    """Protect against infinitesimal heights and associated high velocities
661    """
662
663    #Shortcuts
664    wc = domain.quantities['stage'].centroid_values
665    zc = domain.quantities['elevation'].centroid_values
666    xmomc = domain.quantities['xmomentum'].centroid_values
667    ymomc = domain.quantities['ymomentum'].centroid_values
668    hc = wc - zc  #Water depths at centroids
669
670    #Update
671    #FIXME: Modify accroditg to c-version - or discard altogether.
672    for k in range(domain.number_of_elements):
673
674        if hc[k] < domain.minimum_allowed_height:
675            #Control stage
676            if hc[k] < domain.epsilon:
677                wc[k] = zc[k] # Contain 'lost mass' error
678
679            #Control momentum
680            xmomc[k] = ymomc[k] = 0.0
681
682
683def protect_against_infinitesimal_and_negative_heights_c(domain):
684    """Protect against infinitesimal heights and associated high velocities
685    """
686
687    #Shortcuts
688    wc = domain.quantities['stage'].centroid_values
689    zc = domain.quantities['elevation'].centroid_values
690    xmomc = domain.quantities['xmomentum'].centroid_values
691    ymomc = domain.quantities['ymomentum'].centroid_values
692
693    from shallow_water_ext import protect
694
695    protect(domain.minimum_allowed_height, domain.maximum_allowed_speed, 
696            domain.epsilon, wc, zc, xmomc, ymomc)
697
698
699
700def h_limiter(domain):
701    """Limit slopes for each volume to eliminate artificial variance
702    introduced by e.g. second order extrapolator
703
704    limit on h = w-z
705
706    This limiter depends on two quantities (w,z) so it resides within
707    this module rather than within quantity.py
708    """
709
710    from Numeric import zeros, Float
711
712    N = domain.number_of_elements
713    beta_h = domain.beta_h
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    hvbar = zeros(hv.shape, Float) #h-limited values
725
726    #Find min and max of this and neighbour's centroid values
727    hmax = zeros(hc.shape, Float)
728    hmin = zeros(hc.shape, Float)
729
730    for k in range(N):
731        hmax[k] = hmin[k] = hc[k]
732        for i in range(3):
733            n = domain.neighbours[k,i]
734            if n >= 0:
735                hn = hc[n] #Neighbour's centroid value
736
737                hmin[k] = min(hmin[k], hn)
738                hmax[k] = max(hmax[k], hn)
739
740
741    #Diffences between centroids and maxima/minima
742    dhmax = hmax - hc
743    dhmin = hmin - hc
744
745    #Deltas between vertex and centroid values
746    dh = zeros(hv.shape, Float)
747    for i in range(3):
748        dh[:,i] = hv[:,i] - hc
749
750    #Phi limiter
751    for k in range(N):
752
753        #Find the gradient limiter (phi) across vertices
754        phi = 1.0
755        for i in range(3):
756            r = 1.0
757            if (dh[k,i] > 0): r = dhmax[k]/dh[k,i]
758            if (dh[k,i] < 0): r = dhmin[k]/dh[k,i]
759
760            phi = min( min(r*beta_h, 1), phi )
761
762        #Then update using phi limiter
763        for i in range(3):
764            hvbar[k,i] = hc[k] + phi*dh[k,i]
765
766    return hvbar
767
768
769
770def h_limiter_c(domain):
771    """Limit slopes for each volume to eliminate artificial variance
772    introduced by e.g. second order extrapolator
773
774    limit on h = w-z
775
776    This limiter depends on two quantities (w,z) so it resides within
777    this module rather than within quantity.py
778
779    Wrapper for c-extension
780    """
781
782    from Numeric import zeros, Float
783
784    N = domain.number_of_elements
785    beta_h = domain.beta_h
786
787    #Shortcuts
788    wc = domain.quantities['stage'].centroid_values
789    zc = domain.quantities['elevation'].centroid_values
790    hc = wc - zc
791
792    wv = domain.quantities['stage'].vertex_values
793    zv = domain.quantities['elevation'].vertex_values
794    hv = wv - zv
795
796    #Call C-extension
797    from shallow_water_ext import h_limiter_sw as h_limiter
798    hvbar = h_limiter(domain, hc, hv)
799
800    return hvbar
801
802
803def balance_deep_and_shallow(domain):
804    """Compute linear combination between stage as computed by
805    gradient-limiters limiting using w, and stage computed by
806    gradient-limiters limiting using h (h-limiter).
807    The former takes precedence when heights are large compared to the
808    bed slope while the latter takes precedence when heights are
809    relatively small.  Anything in between is computed as a balanced
810    linear combination in order to avoid numerical disturbances which
811    would otherwise appear as a result of hard switching between
812    modes.
813
814    The h-limiter is always applied irrespective of the order.
815    """
816
817    #Shortcuts
818    wc = domain.quantities['stage'].centroid_values
819    zc = domain.quantities['elevation'].centroid_values
820    hc = wc - zc
821
822    wv = domain.quantities['stage'].vertex_values
823    zv = domain.quantities['elevation'].vertex_values
824    hv = wv-zv
825
826    #Limit h
827    hvbar = h_limiter(domain)
828
829    for k in range(domain.number_of_elements):
830        #Compute maximal variation in bed elevation
831        #  This quantitiy is
832        #    dz = max_i abs(z_i - z_c)
833        #  and it is independent of dimension
834        #  In the 1d case zc = (z0+z1)/2
835        #  In the 2d case zc = (z0+z1+z2)/3
836
837        dz = max(abs(zv[k,0]-zc[k]),
838                 abs(zv[k,1]-zc[k]),
839                 abs(zv[k,2]-zc[k]))
840
841
842        hmin = min( hv[k,:] )
843
844        #Create alpha in [0,1], where alpha==0 means using the h-limited
845        #stage and alpha==1 means using the w-limited stage as
846        #computed by the gradient limiter (both 1st or 2nd order)
847
848        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
849        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
850
851        if dz > 0.0:
852            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
853        else:
854            #Flat bed
855            alpha = 1.0
856
857        #Let
858        #
859        #  wvi be the w-limited stage (wvi = zvi + hvi)
860        #  wvi- be the h-limited state (wvi- = zvi + hvi-)
861        #
862        #
863        #where i=0,1,2 denotes the vertex ids
864        #
865        #Weighted balance between w-limited and h-limited stage is
866        #
867        #  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
868        #
869        #It follows that the updated wvi is
870        #  wvi := zvi + (1-alpha)*hvi- + alpha*hvi
871        #
872        # Momentum is balanced between constant and limited
873
874
875        #for i in range(3):
876        #    wv[k,i] = zv[k,i] + hvbar[k,i]
877
878        #return
879
880        if alpha < 1:
881
882            for i in range(3):
883                wv[k,i] = zv[k,i] + (1-alpha)*hvbar[k,i] + alpha*hv[k,i]
884
885            #Momentums at centroids
886            xmomc = domain.quantities['xmomentum'].centroid_values
887            ymomc = domain.quantities['ymomentum'].centroid_values
888
889            #Momentums at vertices
890            xmomv = domain.quantities['xmomentum'].vertex_values
891            ymomv = domain.quantities['ymomentum'].vertex_values
892
893            # Update momentum as a linear combination of
894            # xmomc and ymomc (shallow) and momentum
895            # from extrapolator xmomv and ymomv (deep).
896            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]
897            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
898
899
900def balance_deep_and_shallow_c(domain):
901    """Wrapper for C implementation
902    """
903
904    #Shortcuts
905    wc = domain.quantities['stage'].centroid_values
906    zc = domain.quantities['elevation'].centroid_values
907    hc = wc - zc
908
909    wv = domain.quantities['stage'].vertex_values
910    zv = domain.quantities['elevation'].vertex_values
911    hv = wv - zv
912
913    #Momentums at centroids
914    xmomc = domain.quantities['xmomentum'].centroid_values
915    ymomc = domain.quantities['ymomentum'].centroid_values
916
917    #Momentums at vertices
918    xmomv = domain.quantities['xmomentum'].vertex_values
919    ymomv = domain.quantities['ymomentum'].vertex_values
920
921    #Limit h
922    hvbar = h_limiter(domain)
923
924    #This is how one would make a first order h_limited value
925    #as in the old balancer (pre 17 Feb 2005):
926    #from Numeric import zeros, Float
927    #hvbar = zeros( (len(hc), 3), Float)
928    #for i in range(3):
929    #    hvbar[:,i] = hc[:]
930
931    from shallow_water_ext import balance_deep_and_shallow
932    balance_deep_and_shallow(wc, zc, hc, wv, zv, hv, hvbar,
933                             xmomc, ymomc, xmomv, ymomv)
934
935
936
937
938###############################################
939#Boundaries - specific to the shallow water wave equation
940class Reflective_boundary(Boundary):
941    """Reflective boundary returns same conserved quantities as
942    those present in its neighbour volume but reflected.
943
944    This class is specific to the shallow water equation as it
945    works with the momentum quantities assumed to be the second
946    and third conserved quantities.
947    """
948
949    def __init__(self, domain = None):
950        Boundary.__init__(self)
951
952        if domain is None:
953            msg = 'Domain must be specified for reflective boundary'
954            raise msg
955
956        #Handy shorthands
957        self.stage   = domain.quantities['stage'].edge_values
958        self.xmom    = domain.quantities['xmomentum'].edge_values
959        self.ymom    = domain.quantities['ymomentum'].edge_values
960        self.normals = domain.normals
961
962        from Numeric import zeros, Float
963        self.conserved_quantities = zeros(3, Float)
964
965    def __repr__(self):
966        return 'Reflective_boundary'
967
968
969    def evaluate(self, vol_id, edge_id):
970        """Reflective boundaries reverses the outward momentum
971        of the volume they serve.
972        """
973
974        q = self.conserved_quantities
975        q[0] = self.stage[vol_id, edge_id]
976        q[1] = self.xmom[vol_id, edge_id]
977        q[2] = self.ymom[vol_id, edge_id]
978
979        normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
980
981
982        r = rotate(q, normal, direction = 1)
983        r[1] = -r[1]
984        q = rotate(r, normal, direction = -1)
985
986        return q
987
988
989
990class Transmissive_Momentum_Set_Stage_boundary(Boundary):
991    """Returns same momentum conserved quantities as
992    those present in its neighbour volume. Sets stage
993
994    Underlying domain must be specified when boundary is instantiated
995    """
996
997    def __init__(self, domain = None, function=None):
998        Boundary.__init__(self)
999
1000        if domain is None:
1001            msg = 'Domain must be specified for this type boundary'
1002            raise msg
1003
1004        if function is None:
1005            msg = 'Function must be specified for this type boundary'
1006            raise msg
1007
1008        self.domain   = domain
1009        self.function = function
1010
1011    def __repr__(self):
1012        return 'Transmissive_Momentum_Set_Stage_boundary(%s)' %self.domain
1013
1014    def evaluate(self, vol_id, edge_id):
1015        """Transmissive Momentum Set Stage boundaries return the edge momentum
1016        values of the volume they serve.
1017        """
1018
1019        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
1020        value = self.function(self.domain.time)
1021        q[0] = value[0]
1022        return q
1023
1024
1025        #FIXME: Consider this (taken from File_boundary) to allow
1026        #spatial variation
1027        #if vol_id is not None and edge_id is not None:
1028        #    i = self.boundary_indices[ vol_id, edge_id ]
1029        #    return self.F(t, point_id = i)
1030        #else:
1031        #    return self.F(t)
1032
1033
1034
1035class Dirichlet_Discharge_boundary(Boundary):
1036    """
1037    Sets stage (stage0)
1038    Sets momentum (wh0) in the inward normal direction.
1039
1040    Underlying domain must be specified when boundary is instantiated
1041    """
1042
1043    def __init__(self, domain = None, stage0=None, wh0=None):
1044        Boundary.__init__(self)
1045
1046        if domain is None:
1047            msg = 'Domain must be specified for this type boundary'
1048            raise msg
1049
1050        if stage0 is None:
1051            raise 'set stage'
1052
1053        if wh0 is None:
1054            wh0 = 0.0
1055
1056        self.domain   = domain
1057        self.stage0  = stage0
1058        self.wh0 = wh0
1059
1060    def __repr__(self):
1061        return 'Dirichlet_Discharge_boundary(%s)' %self.domain
1062
1063    def evaluate(self, vol_id, edge_id):
1064        """Set discharge in the (inward) normal direction
1065        """
1066
1067        normal = self.domain.get_normal(vol_id,edge_id)
1068        q = [self.stage0, -self.wh0*normal[0], -self.wh0*normal[1]]
1069        return q
1070
1071
1072        #FIXME: Consider this (taken from File_boundary) to allow
1073        #spatial variation
1074        #if vol_id is not None and edge_id is not None:
1075        #    i = self.boundary_indices[ vol_id, edge_id ]
1076        #    return self.F(t, point_id = i)
1077        #else:
1078        #    return self.F(t)
1079
1080
1081
1082#class Spatio_temporal_boundary(Boundary):
1083#    """The spatio-temporal boundary, reads values for the conserved
1084#    quantities from an sww NetCDF file, and returns interpolated values
1085#    at the midpoints of each associated boundaty segment.
1086#    Time dependency is interpolated linearly as in util.File_function.#
1087#
1088#    Example:
1089#    Bf = Spatio_temporal_boundary('source_file.sww', domain)
1090#
1091#    """
1092Spatio_temporal_boundary = File_boundary
1093
1094
1095
1096
1097#########################
1098#Standard forcing terms:
1099#
1100def gravity(domain):
1101    """Apply gravitational pull in the presence of bed slope
1102    """
1103
1104    from Numeric import zeros, Float, array, sum
1105
1106    xmom = domain.quantities['xmomentum'].explicit_update
1107    ymom = domain.quantities['ymomentum'].explicit_update
1108
1109    Stage = domain.quantities['stage']
1110    Elevation = domain.quantities['elevation']
1111    h = Stage.edge_values - Elevation.edge_values
1112    v = Elevation.vertex_values
1113
1114    x = domain.get_vertex_coordinates()
1115    g = domain.g
1116
1117    for k in range(domain.number_of_elements):
1118        avg_h = sum( h[k,:] )/3
1119
1120        #Compute bed slope
1121        x0, y0, x1, y1, x2, y2 = x[k,:]
1122        z0, z1, z2 = v[k,:]
1123
1124        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
1125
1126        #Update momentum
1127        xmom[k] += -g*zx*avg_h
1128        ymom[k] += -g*zy*avg_h
1129
1130
1131def gravity_c(domain):
1132    """Wrapper calling C version
1133    """
1134
1135    xmom = domain.quantities['xmomentum'].explicit_update
1136    ymom = domain.quantities['ymomentum'].explicit_update
1137
1138    Stage = domain.quantities['stage']
1139    Elevation = domain.quantities['elevation']
1140    h = Stage.edge_values - Elevation.edge_values
1141    v = Elevation.vertex_values
1142
1143    x = domain.get_vertex_coordinates()
1144    g = domain.g
1145
1146
1147    from shallow_water_ext import gravity
1148    gravity(g, h, v, x, xmom, ymom)
1149
1150
1151def manning_friction(domain):
1152    """Apply (Manning) friction to water momentum
1153    """
1154
1155    from math import sqrt
1156
1157    w = domain.quantities['stage'].centroid_values
1158    z = domain.quantities['elevation'].centroid_values
1159    h = w-z
1160
1161    uh = domain.quantities['xmomentum'].centroid_values
1162    vh = domain.quantities['ymomentum'].centroid_values
1163    eta = domain.quantities['friction'].centroid_values
1164
1165    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1166    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1167
1168    N = domain.number_of_elements
1169    eps = domain.minimum_allowed_height
1170    g = domain.g
1171
1172    for k in range(N):
1173        if eta[k] >= eps:
1174            if h[k] >= eps:
1175                S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
1176                S /= h[k]**(7.0/3)
1177
1178                #Update momentum
1179                xmom_update[k] += S*uh[k]
1180                ymom_update[k] += S*vh[k]
1181
1182
1183def manning_friction_c(domain):
1184    """Wrapper for c version
1185    """
1186
1187
1188    xmom = domain.quantities['xmomentum']
1189    ymom = domain.quantities['ymomentum']
1190
1191    w = domain.quantities['stage'].centroid_values
1192    z = domain.quantities['elevation'].centroid_values
1193
1194    uh = xmom.centroid_values
1195    vh = ymom.centroid_values
1196    eta = domain.quantities['friction'].centroid_values
1197
1198    xmom_update = xmom.semi_implicit_update
1199    ymom_update = ymom.semi_implicit_update
1200
1201    N = domain.number_of_elements
1202    eps = domain.minimum_allowed_height
1203    g = domain.g
1204
1205    from shallow_water_ext import manning_friction
1206    manning_friction(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
1207
1208
1209def linear_friction(domain):
1210    """Apply linear friction to water momentum
1211
1212    Assumes quantity: 'linear_friction' to be present
1213    """
1214
1215    from math import sqrt
1216
1217    w = domain.quantities['stage'].centroid_values
1218    z = domain.quantities['elevation'].centroid_values
1219    h = w-z
1220
1221    uh = domain.quantities['xmomentum'].centroid_values
1222    vh = domain.quantities['ymomentum'].centroid_values
1223    tau = domain.quantities['linear_friction'].centroid_values
1224
1225    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1226    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1227
1228    N = domain.number_of_elements
1229    eps = domain.minimum_allowed_height
1230    g = domain.g #Not necessary? Why was this added?
1231
1232    for k in range(N):
1233        if tau[k] >= eps:
1234            if h[k] >= eps:
1235                S = -tau[k]/h[k]
1236
1237                #Update momentum
1238                xmom_update[k] += S*uh[k]
1239                ymom_update[k] += S*vh[k]
1240
1241
1242
1243def check_forcefield(f):
1244    """Check that f is either
1245    1: a callable object f(t,x,y), where x and y are vectors
1246       and that it returns an array or a list of same length
1247       as x and y
1248    2: a scalar
1249    """
1250
1251    from Numeric import ones, Float, array
1252
1253
1254    if callable(f):
1255        N = 3
1256        x = ones(3, Float)
1257        y = ones(3, Float)
1258        try:
1259            q = f(1.0, x=x, y=y)
1260        except Exception, e:
1261            msg = 'Function %s could not be executed:\n%s' %(f, e)
1262            #FIXME: Reconsider this semantics
1263            raise msg
1264
1265        try:
1266            q = array(q).astype(Float)
1267        except:
1268            msg = 'Return value from vector function %s could ' %f
1269            msg += 'not be converted into a Numeric array of floats.\n'
1270            msg += 'Specified function should return either list or array.'
1271            raise msg
1272
1273        #Is this really what we want?
1274        msg = 'Return vector from function %s ' %f
1275        msg += 'must have same lenght as input vectors'
1276        assert len(q) == N, msg
1277
1278    else:
1279        try:
1280            f = float(f)
1281        except:
1282            msg = 'Force field %s must be either a scalar' %f
1283            msg += ' or a vector function'
1284            raise msg
1285    return f
1286
1287
1288class Wind_stress:
1289    """Apply wind stress to water momentum in terms of
1290    wind speed [m/s] and wind direction [degrees]
1291    """
1292
1293    def __init__(self, *args, **kwargs):
1294        """Initialise windfield from wind speed s [m/s]
1295        and wind direction phi [degrees]
1296
1297        Inputs v and phi can be either scalars or Python functions, e.g.
1298
1299        W = Wind_stress(10, 178)
1300
1301        #FIXME - 'normal' degrees are assumed for now, i.e. the
1302        vector (1,0) has zero degrees.
1303        We may need to convert from 'compass' degrees later on and also
1304        map from True north to grid north.
1305
1306        Arguments can also be Python functions of t,x,y as in
1307
1308        def speed(t,x,y):
1309            ...
1310            return s
1311
1312        def angle(t,x,y):
1313            ...
1314            return phi
1315
1316        where x and y are vectors.
1317
1318        and then pass the functions in
1319
1320        W = Wind_stress(speed, angle)
1321
1322        The instantiated object W can be appended to the list of
1323        forcing_terms as in
1324
1325        Alternatively, one vector valued function for (speed, angle)
1326        can be applied, providing both quantities simultaneously.
1327        As in
1328        W = Wind_stress(F), where returns (speed, angle) for each t.
1329
1330        domain.forcing_terms.append(W)
1331        """
1332
1333        from config import rho_a, rho_w, eta_w
1334        from Numeric import array, Float
1335
1336        if len(args) == 2:
1337            s = args[0]
1338            phi = args[1]
1339        elif len(args) == 1:
1340            #Assume vector function returning (s, phi)(t,x,y)
1341            vector_function = args[0]
1342            s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
1343            phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
1344        else:
1345           #Assume info is in 2 keyword arguments
1346
1347           if len(kwargs) == 2:
1348               s = kwargs['s']
1349               phi = kwargs['phi']
1350           else:
1351               raise 'Assumes two keyword arguments: s=..., phi=....'
1352
1353        self.speed = check_forcefield(s)
1354        self.phi = check_forcefield(phi)
1355
1356        self.const = eta_w*rho_a/rho_w
1357
1358
1359    def __call__(self, domain):
1360        """Evaluate windfield based on values found in domain
1361        """
1362
1363        from math import pi, cos, sin, sqrt
1364        from Numeric import Float, ones, ArrayType
1365
1366        xmom_update = domain.quantities['xmomentum'].explicit_update
1367        ymom_update = domain.quantities['ymomentum'].explicit_update
1368
1369        N = domain.number_of_elements
1370        t = domain.time
1371
1372        if callable(self.speed):
1373            xc = domain.get_centroid_coordinates()
1374            s_vec = self.speed(t, xc[:,0], xc[:,1])
1375        else:
1376            #Assume s is a scalar
1377
1378            try:
1379                s_vec = self.speed * ones(N, Float)
1380            except:
1381                msg = 'Speed must be either callable or a scalar: %s' %self.s
1382                raise msg
1383
1384
1385        if callable(self.phi):
1386            xc = domain.get_centroid_coordinates()
1387            phi_vec = self.phi(t, xc[:,0], xc[:,1])
1388        else:
1389            #Assume phi is a scalar
1390
1391            try:
1392                phi_vec = self.phi * ones(N, Float)
1393            except:
1394                msg = 'Angle must be either callable or a scalar: %s' %self.phi
1395                raise msg
1396
1397        assign_windfield_values(xmom_update, ymom_update,
1398                                s_vec, phi_vec, self.const)
1399
1400
1401def assign_windfield_values(xmom_update, ymom_update,
1402                            s_vec, phi_vec, const):
1403    """Python version of assigning wind field to update vectors.
1404    A c version also exists (for speed)
1405    """
1406    from math import pi, cos, sin, sqrt
1407
1408    N = len(s_vec)
1409    for k in range(N):
1410        s = s_vec[k]
1411        phi = phi_vec[k]
1412
1413        #Convert to radians
1414        phi = phi*pi/180
1415
1416        #Compute velocity vector (u, v)
1417        u = s*cos(phi)
1418        v = s*sin(phi)
1419
1420        #Compute wind stress
1421        S = const * sqrt(u**2 + v**2)
1422        xmom_update[k] += S*u
1423        ymom_update[k] += S*v
1424
1425
1426
1427##############################
1428#OBSOLETE STUFF
1429
1430def balance_deep_and_shallow_old(domain):
1431    """Compute linear combination between stage as computed by
1432    gradient-limiters and stage computed as constant height above bed.
1433    The former takes precedence when heights are large compared to the
1434    bed slope while the latter takes precedence when heights are
1435    relatively small.  Anything in between is computed as a balanced
1436    linear combination in order to avoid numerical disturbances which
1437    would otherwise appear as a result of hard switching between
1438    modes.
1439    """
1440
1441    #OBSOLETE
1442
1443    #Shortcuts
1444    wc = domain.quantities['stage'].centroid_values
1445    zc = domain.quantities['elevation'].centroid_values
1446    hc = wc - zc
1447
1448    wv = domain.quantities['stage'].vertex_values
1449    zv = domain.quantities['elevation'].vertex_values
1450    hv = wv-zv
1451
1452
1453    #Computed linear combination between constant stages and and
1454    #stages parallel to the bed elevation.
1455    for k in range(domain.number_of_elements):
1456        #Compute maximal variation in bed elevation
1457        #  This quantitiy is
1458        #    dz = max_i abs(z_i - z_c)
1459        #  and it is independent of dimension
1460        #  In the 1d case zc = (z0+z1)/2
1461        #  In the 2d case zc = (z0+z1+z2)/3
1462
1463        dz = max(abs(zv[k,0]-zc[k]),
1464                 abs(zv[k,1]-zc[k]),
1465                 abs(zv[k,2]-zc[k]))
1466
1467
1468        hmin = min( hv[k,:] )
1469
1470        #Create alpha in [0,1], where alpha==0 means using shallow
1471        #first order scheme and alpha==1 means using the stage w as
1472        #computed by the gradient limiter (1st or 2nd order)
1473        #
1474        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
1475        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
1476
1477        if dz > 0.0:
1478            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
1479        else:
1480            #Flat bed
1481            alpha = 1.0
1482
1483        #Weighted balance between stage parallel to bed elevation
1484        #(wvi = zvi + hc) and stage as computed by 1st or 2nd
1485        #order gradient limiter
1486        #(wvi = zvi + hvi) where i=0,1,2 denotes the vertex ids
1487        #
1488        #It follows that the updated wvi is
1489        #  wvi := (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) =
1490        #  zvi + hc + alpha*(hvi - hc)
1491        #
1492        #Note that hvi = zc+hc-zvi in the first order case (constant).
1493
1494        if alpha < 1:
1495            for i in range(3):
1496                wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k])
1497
1498
1499            #Momentums at centroids
1500            xmomc = domain.quantities['xmomentum'].centroid_values
1501            ymomc = domain.quantities['ymomentum'].centroid_values
1502
1503            #Momentums at vertices
1504            xmomv = domain.quantities['xmomentum'].vertex_values
1505            ymomv = domain.quantities['ymomentum'].vertex_values
1506
1507            # Update momentum as a linear combination of
1508            # xmomc and ymomc (shallow) and momentum
1509            # from extrapolator xmomv and ymomv (deep).
1510            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]
1511            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
1512
1513
1514
1515
1516
1517###########################
1518###########################
1519#Geometries
1520
1521
1522#FIXME: Rethink this way of creating values.
1523
1524
1525class Weir:
1526    """Set a bathymetry for weir with a hole and a downstream gutter
1527    x,y are assumed to be in the unit square
1528    """
1529
1530    def __init__(self, stage):
1531        self.inflow_stage = stage
1532
1533    def __call__(self, x, y):
1534        from Numeric import zeros, Float
1535        from math import sqrt
1536
1537        N = len(x)
1538        assert N == len(y)
1539
1540        z = zeros(N, Float)
1541        for i in range(N):
1542            z[i] = -x[i]/2  #General slope
1543
1544            #Flattish bit to the left
1545            if x[i] < 0.3:
1546                z[i] = -x[i]/10
1547
1548            #Weir
1549            if x[i] >= 0.3 and x[i] < 0.4:
1550                z[i] = -x[i]+0.9
1551
1552            #Dip
1553            x0 = 0.6
1554            #depth = -1.3
1555            depth = -1.0
1556            #plateaux = -0.9
1557            plateaux = -0.6
1558            if y[i] < 0.7:
1559                if x[i] > x0 and x[i] < 0.9:
1560                    z[i] = depth
1561
1562                #RHS plateaux
1563                if x[i] >= 0.9:
1564                    z[i] = plateaux
1565
1566
1567            elif y[i] >= 0.7 and y[i] < 1.5:
1568                #Restrict and deepen
1569                if x[i] >= x0 and x[i] < 0.8:
1570                    z[i] = depth-(y[i]/3-0.3)
1571                    #z[i] = depth-y[i]/5
1572                    #z[i] = depth
1573                elif x[i] >= 0.8:
1574                    #RHS plateaux
1575                    z[i] = plateaux
1576
1577            elif y[i] >= 1.5:
1578                if x[i] >= x0 and x[i] < 0.8 + (y[i]-1.5)/1.2:
1579                    #Widen up and stay at constant depth
1580                    z[i] = depth-1.5/5
1581                elif x[i] >= 0.8 + (y[i]-1.5)/1.2:
1582                    #RHS plateaux
1583                    z[i] = plateaux
1584
1585
1586            #Hole in weir (slightly higher than inflow condition)
1587            if x[i] >= 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
1588                z[i] = -x[i]+self.inflow_stage + 0.02
1589
1590            #Channel behind weir
1591            x0 = 0.5
1592            if x[i] >= 0.4 and x[i] < x0 and y[i] > 0.2 and y[i] < 0.4:
1593                z[i] = -x[i]+self.inflow_stage + 0.02
1594
1595            if x[i] >= x0 and x[i] < 0.6 and y[i] > 0.2 and y[i] < 0.4:
1596                #Flatten it out towards the end
1597                z[i] = -x0+self.inflow_stage + 0.02 + (x0-x[i])/5
1598
1599            #Hole to the east
1600            x0 = 1.1; y0 = 0.35
1601            #if x[i] < -0.2 and y < 0.5:
1602            if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
1603                z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-1.0
1604
1605            #Tiny channel draining hole
1606            if x[i] >= 1.14 and x[i] < 1.2 and y[i] >= 0.4 and y[i] < 0.6:
1607                z[i] = -0.9 #North south
1608
1609            if x[i] >= 0.9 and x[i] < 1.18 and y[i] >= 0.58 and y[i] < 0.65:
1610                z[i] = -1.0 + (x[i]-0.9)/3 #East west
1611
1612
1613
1614            #Stuff not in use
1615
1616            #Upward slope at inlet to the north west
1617            #if x[i] < 0.0: # and y[i] > 0.5:
1618            #    #z[i] = -y[i]+0.5  #-x[i]/2
1619            #    z[i] = x[i]/4 - y[i]**2 + 0.5
1620
1621            #Hole to the west
1622            #x0 = -0.4; y0 = 0.35 # center
1623            #if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
1624            #    z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-0.2
1625
1626
1627
1628
1629
1630        return z/2
1631
1632class Weir_simple:
1633    """Set a bathymetry for weir with a hole and a downstream gutter
1634    x,y are assumed to be in the unit square
1635    """
1636
1637    def __init__(self, stage):
1638        self.inflow_stage = stage
1639
1640    def __call__(self, x, y):
1641        from Numeric import zeros, Float
1642
1643        N = len(x)
1644        assert N == len(y)
1645
1646        z = zeros(N, Float)
1647        for i in range(N):
1648            z[i] = -x[i]  #General slope
1649
1650            #Flat bit to the left
1651            if x[i] < 0.3:
1652                z[i] = -x[i]/10  #General slope
1653
1654            #Weir
1655            if x[i] > 0.3 and x[i] < 0.4:
1656                z[i] = -x[i]+0.9
1657
1658            #Dip
1659            if x[i] > 0.6 and x[i] < 0.9:
1660                z[i] = -x[i]-0.5  #-y[i]/5
1661
1662            #Hole in weir (slightly higher than inflow condition)
1663            if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
1664                z[i] = -x[i]+self.inflow_stage + 0.05
1665
1666
1667        return z/2
1668
1669
1670
1671class Constant_stage:
1672    """Set an initial condition with constant stage
1673    """
1674    def __init__(self, s):
1675        self.s = s
1676
1677    def __call__(self, x, y):
1678        return self.s
1679
1680class Constant_height:
1681    """Set an initial condition with constant water height, e.g
1682    stage s = z+h
1683    """
1684
1685    def __init__(self, W, h):
1686        self.W = W
1687        self.h = h
1688
1689    def __call__(self, x, y):
1690        if self.W is None:
1691            from Numeric import ones, Float
1692            return self.h*ones(len(x), Float)
1693        else:
1694            return self.W(x,y) + self.h
1695
1696
1697
1698
1699def compute_fluxes_python(domain):
1700    """Compute all fluxes and the timestep suitable for all volumes
1701    in domain.
1702
1703    Compute total flux for each conserved quantity using "flux_function"
1704
1705    Fluxes across each edge are scaled by edgelengths and summed up
1706    Resulting flux is then scaled by area and stored in
1707    explicit_update for each of the three conserved quantities
1708    stage, xmomentum and ymomentum
1709
1710    The maximal allowable speed computed by the flux_function for each volume
1711    is converted to a timestep that must not be exceeded. The minimum of
1712    those is computed as the next overall timestep.
1713
1714    Post conditions:
1715      domain.explicit_update is reset to computed flux values
1716      domain.timestep is set to the largest step satisfying all volumes.
1717    """
1718
1719    import sys
1720    from Numeric import zeros, Float
1721
1722    N = domain.number_of_elements
1723
1724    #Shortcuts
1725    Stage = domain.quantities['stage']
1726    Xmom = domain.quantities['xmomentum']
1727    Ymom = domain.quantities['ymomentum']
1728    Bed = domain.quantities['elevation']
1729
1730    #Arrays
1731    stage = Stage.edge_values
1732    xmom =  Xmom.edge_values
1733    ymom =  Ymom.edge_values
1734    bed =   Bed.edge_values
1735
1736    stage_bdry = Stage.boundary_values
1737    xmom_bdry =  Xmom.boundary_values
1738    ymom_bdry =  Ymom.boundary_values
1739
1740    flux = zeros((N,3), Float) #Work array for summing up fluxes
1741
1742    #Loop
1743    timestep = float(sys.maxint)
1744    for k in range(N):
1745
1746        for i in range(3):
1747            #Quantities inside volume facing neighbour i
1748            ql = [stage[k, i], xmom[k, i], ymom[k, i]]
1749            zl = bed[k, i]
1750
1751            #Quantities at neighbour on nearest face
1752            n = domain.neighbours[k,i]
1753            if n < 0:
1754                m = -n-1 #Convert negative flag to index
1755                qr = [stage_bdry[m], xmom_bdry[m], ymom_bdry[m]]
1756                zr = zl #Extend bed elevation to boundary
1757            else:
1758                m = domain.neighbour_edges[k,i]
1759                qr = [stage[n, m], xmom[n, m], ymom[n, m]]
1760                zr = bed[n, m]
1761
1762
1763            #Outward pointing normal vector
1764            normal = domain.normals[k, 2*i:2*i+2]
1765
1766            #Flux computation using provided function
1767            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
1768
1769            flux[k,:] = edgeflux
1770
1771    return flux
1772
1773
1774
1775
1776
1777
1778
1779##############################################
1780#Initialise module
1781
1782
1783from utilities import compile
1784if compile.can_use_C_extension('shallow_water_ext.c'):
1785    #Replace python version with c implementations
1786
1787    from shallow_water_ext import rotate, assign_windfield_values
1788    compute_fluxes = compute_fluxes_c
1789    extrapolate_second_order_sw=extrapolate_second_order_sw_c
1790    gravity = gravity_c
1791    manning_friction = manning_friction_c
1792    h_limiter = h_limiter_c
1793    balance_deep_and_shallow = balance_deep_and_shallow_c
1794    protect_against_infinitesimal_and_negative_heights = protect_against_infinitesimal_and_negative_heights_c
1795
1796
1797    #distribute_to_vertices_and_edges = distribute_to_vertices_and_edges_c #(like MH's)
1798
1799
1800
1801#Optimisation with psyco
1802from config import use_psyco
1803if use_psyco:
1804    try:
1805        import psyco
1806    except:
1807        import os
1808        if os.name == 'posix' and os.uname()[4] == 'x86_64':
1809            pass
1810            #Psyco isn't supported on 64 bit systems, but it doesn't matter
1811        else:
1812            msg = 'WARNING: psyco (speedup) could not import'+\
1813                  ', you may want to consider installing it'
1814            print msg
1815    else:
1816        psyco.bind(Domain.distribute_to_vertices_and_edges)
1817        psyco.bind(Domain.compute_fluxes)
1818
1819if __name__ == "__main__":
1820    pass
1821
1822# Profiling stuff
1823import profile
1824profiler = profile.Profile()
Note: See TracBrowser for help on using the repository browser.