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

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

Refactored limit2007 to tight_slope_limiters (using eclipse)

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