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

Last change on this file since 1697 was 1697, checked in by steve, 19 years ago

Found problem with File_Boundary as used in validation test LWRU2. Have setup new BC Transmissive_Momentum_Set _Stage

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