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

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