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

Last change on this file since 3876 was 3876, checked in by steve, 17 years ago

Added parameter alpha_balance to config.py to deal with large jumps in bed

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