source: inundation/pyvolution/shallow_water.py @ 1995

Last change on this file since 1995 was 1922, checked in by duncan, 19 years ago

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