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

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

Implemented dry-cell exclusion for the linear reconstruction step.
Time for the okushiri performance example reduced the running time (on my GA windows box)
of that step from 15.3s to 11.6s - an almost 25% improvement.
With a total running time of this example of about 33s before the optimisation this translates to just over 10% total improvement of the running time of the evolve loop.

With more dry land, this optimisation will be more important.

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