source: anuga_core/source/anuga/pyvolution/shallow_water.py @ 3532

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

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