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

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

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