source: inundation/ga/storm_surge/pyvolution/shallow_water.py @ 1581

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