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

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

Added discoverability keywords into docstring

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