source: inundation-numpy-branch/pyvolution/shallow_water.py @ 3514

Last change on this file since 3514 was 3514, checked in by duncan, 18 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

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