source: anuga_core/source/anuga/shallow_water/shallow_water_domain.py @ 3689

Last change on this file since 3689 was 3689, checked in by ole, 17 years ago

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