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

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

Working on 2nd order time stepping

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