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

Last change on this file since 6647 was 6647, checked in by ted, 15 years ago

Started function for volumetric test

  • Property svn:keywords set to LastChangedDate LastChangedRevision LastChangedBy HeadURL Id
File size: 75.1 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
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: 6647 $)
70ModifiedBy:
71    $Author: ted $
72    $Date: 2009-03-27 04:32:36 +0000 (Fri, 27 Mar 2009) $
73
74"""
75
76# Subversion keywords:
77#
78# $LastChangedDate: 2009-03-27 04:32:36 +0000 (Fri, 27 Mar 2009) $
79# $LastChangedRevision: 6647 $
80# $LastChangedBy: ted $
81
82import Numeric as num
83
84from anuga.abstract_2d_finite_volumes.neighbour_mesh import segment_midpoints
85from anuga.abstract_2d_finite_volumes.domain import Domain as Generic_Domain
86from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
87     import Boundary
88from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
89     import File_boundary
90from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
91     import Dirichlet_boundary
92from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
93     import Time_boundary
94from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
95     import Transmissive_boundary
96from anuga.pmesh.mesh_interface import create_mesh_from_regions
97from anuga.utilities.numerical_tools import gradient, mean, ensure_numeric
98from anuga.geospatial_data.geospatial_data import ensure_geospatial
99
100from anuga.config import minimum_storable_height
101from anuga.config import minimum_allowed_height, maximum_allowed_speed
102from anuga.config import g, epsilon, beta_w, beta_w_dry,\
103     beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, tight_slope_limiters
104from anuga.config import alpha_balance
105from anuga.config import optimise_dry_cells
106from anuga.config import optimised_gradient_limiter
107from anuga.config import use_edge_limiter
108from anuga.config import use_centroid_velocities
109from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
110
111from anuga.fit_interpolate.interpolate import Modeltime_too_late, Modeltime_too_early
112
113from anuga.utilities.polygon import inside_polygon, polygon_area, is_inside_polygon
114
115from types import IntType, FloatType
116from warnings import warn
117
118
119
120#
121# Shallow water domain
122#---------------------
123class Domain(Generic_Domain):
124
125    conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
126    other_quantities = ['elevation', 'friction']
127   
128    def __init__(self,
129                 coordinates=None,
130                 vertices=None,
131                 boundary=None,
132                 tagged_elements=None,
133                 geo_reference=None,
134                 use_inscribed_circle=False,
135                 mesh_filename=None,
136                 use_cache=False,
137                 verbose=False,
138                 full_send_dict=None,
139                 ghost_recv_dict=None,
140                 processor=0,
141                 numproc=1,
142                 number_of_full_nodes=None,
143                 number_of_full_triangles=None):
144
145
146        other_quantities = ['elevation', 'friction']
147        Generic_Domain.__init__(self,
148                                coordinates,
149                                vertices,
150                                boundary,
151                                Domain.conserved_quantities,
152                                Domain.other_quantities,
153                                tagged_elements,
154                                geo_reference,
155                                use_inscribed_circle,
156                                mesh_filename,
157                                use_cache,
158                                verbose,
159                                full_send_dict,
160                                ghost_recv_dict,
161                                processor,
162                                numproc,
163                                number_of_full_nodes=number_of_full_nodes,
164                                number_of_full_triangles=number_of_full_triangles) 
165
166
167        self.set_minimum_allowed_height(minimum_allowed_height)
168       
169        self.maximum_allowed_speed = maximum_allowed_speed
170        self.g = g
171        self.beta_w      = beta_w
172        self.beta_w_dry  = beta_w_dry
173        self.beta_uh     = beta_uh
174        self.beta_uh_dry = beta_uh_dry
175        self.beta_vh     = beta_vh
176        self.beta_vh_dry = beta_vh_dry
177        self.alpha_balance = alpha_balance
178
179        self.tight_slope_limiters = tight_slope_limiters
180        self.optimise_dry_cells = optimise_dry_cells
181
182        self.forcing_terms.append(manning_friction_implicit)
183        self.forcing_terms.append(gravity)
184
185        # Stored output
186        self.store = True
187        self.format = 'sww'
188        self.set_store_vertices_uniquely(False)
189        self.minimum_storable_height = minimum_storable_height
190        self.quantities_to_be_stored = ['stage','xmomentum','ymomentum']
191
192        # Limiters
193        self.use_edge_limiter = use_edge_limiter
194        self.optimised_gradient_limiter = optimised_gradient_limiter
195        self.use_centroid_velocities = use_centroid_velocities
196
197
198    def set_all_limiters(self, beta):
199        """Shorthand to assign one constant value [0,1[ to all limiters.
200        0 Corresponds to first order, where as larger values make use of
201        the second order scheme.
202        """
203
204        self.beta_w      = beta
205        self.beta_w_dry  = beta
206        self.quantities['stage'].beta = beta
207       
208        self.beta_uh     = beta
209        self.beta_uh_dry = beta
210        self.quantities['xmomentum'].beta = beta
211       
212        self.beta_vh     = beta
213        self.beta_vh_dry = beta
214        self.quantities['ymomentum'].beta = beta
215       
216       
217
218    def set_store_vertices_uniquely(self, flag, reduction=None):
219        """Decide whether vertex values should be stored uniquely as
220        computed in the model or whether they should be reduced to one
221        value per vertex using self.reduction.
222        """
223
224        # FIXME (Ole): how about using the word continuous vertex values?
225        self.smooth = not flag
226
227        # Reduction operation for get_vertex_values
228        if reduction is None:
229            self.reduction = mean
230            #self.reduction = min  #Looks better near steep slopes
231
232
233    def set_minimum_storable_height(self, minimum_storable_height):
234        """
235        Set the minimum depth that will be recognised when writing
236        to an sww file. This is useful for removing thin water layers
237        that seems to be caused by friction creep.
238
239        The minimum allowed sww depth is in meters.
240        """
241        self.minimum_storable_height = minimum_storable_height
242
243
244    def set_minimum_allowed_height(self, minimum_allowed_height):
245        """
246        Set the minimum depth that will be recognised in the numerical
247        scheme
248
249        The minimum allowed depth is in meters.
250
251        The parameter H0 (Minimal height for flux computation)
252        is also set by this function
253        """
254
255        #FIXME (Ole): rename H0 to minimum_allowed_height_in_flux_computation
256
257        #FIXME (Ole): Maybe use histogram to identify isolated extreme speeds
258        #and deal with them adaptively similarly to how we used to use 1 order
259        #steps to recover.
260        self.minimum_allowed_height = minimum_allowed_height
261        self.H0 = minimum_allowed_height   
262       
263
264    def set_maximum_allowed_speed(self, maximum_allowed_speed):
265        """
266        Set the maximum particle speed that is allowed in water
267        shallower than minimum_allowed_height. This is useful for
268        controlling speeds in very thin layers of water and at the same time
269        allow some movement avoiding pooling of water.
270
271        """
272        self.maximum_allowed_speed = maximum_allowed_speed
273
274
275    def set_points_file_block_line_size(self,points_file_block_line_size):
276        """
277        Set the minimum depth that will be recognised when writing
278        to an sww file. This is useful for removing thin water layers
279        that seems to be caused by friction creep.
280
281        The minimum allowed sww depth is in meters.
282        """
283        self.points_file_block_line_size = points_file_block_line_size
284       
285       
286    def set_quantities_to_be_stored(self, q):
287        """Specify which quantities will be stored in the sww file.
288
289        q must be either:
290          - the name of a quantity
291          - a list of quantity names
292          - None
293
294        In the two first cases, the named quantities will be stored at
295        each yieldstep (This is in addition to the quantities elevation
296        and friction)
297       
298        If q is None, storage will be switched off altogether.
299        """
300
301
302        if q is None:
303            self.quantities_to_be_stored = []
304            self.store = False
305            return
306
307        if isinstance(q, basestring):
308            q = [q] # Turn argument into a list
309
310        # Check correcness
311        for quantity_name in q:
312            msg = 'Quantity %s is not a valid conserved quantity'\
313                  %quantity_name
314           
315            assert quantity_name in self.conserved_quantities, msg
316
317        self.quantities_to_be_stored = q
318
319
320
321    def get_wet_elements(self, indices=None):
322        """Return indices for elements where h > minimum_allowed_height
323
324        Optional argument:
325            indices is the set of element ids that the operation applies to.
326
327        Usage:
328            indices = get_wet_elements()
329
330        Note, centroid values are used for this operation           
331        """
332
333        # Water depth below which it is considered to be 0 in the model
334        # FIXME (Ole): Allow this to be specified as a keyword argument as well
335        from anuga.config import minimum_allowed_height
336
337       
338        elevation = self.get_quantity('elevation').\
339                    get_values(location='centroids', indices=indices)
340        stage = self.get_quantity('stage').\
341                get_values(location='centroids', indices=indices)       
342        depth = stage - elevation
343
344        # Select indices for which depth > 0
345        wet_indices = num.compress(depth > minimum_allowed_height,
346                                   num.arange(len(depth)))
347        return wet_indices
348
349
350    def get_maximum_inundation_elevation(self, indices=None):
351        """Return highest elevation where h > 0
352
353        Optional argument:
354            indices is the set of element ids that the operation applies to.
355
356        Usage:
357            q = get_maximum_inundation_elevation()
358
359        Note, centroid values are used for this operation           
360        """
361
362        wet_elements = self.get_wet_elements(indices)
363        return self.get_quantity('elevation').\
364               get_maximum_value(indices=wet_elements)
365
366
367    def get_maximum_inundation_location(self, indices=None):
368        """Return location of highest elevation where h > 0
369
370        Optional argument:
371            indices is the set of element ids that the operation applies to.
372
373        Usage:
374            q = get_maximum_inundation_location()
375
376        Note, centroid values are used for this operation           
377        """
378
379        wet_elements = self.get_wet_elements(indices)
380        return self.get_quantity('elevation').\
381               get_maximum_location(indices=wet_elements)   
382               
383               
384               
385               
386    def get_flow_through_cross_section(self, polyline,
387                                       verbose=False):               
388        """Get the total flow through an arbitrary poly line.       
389       
390        This is a run-time equivalent of the function with same name
391        in data_manager.py
392       
393        Input:
394            polyline: Representation of desired cross section - it may contain
395                      multiple sections allowing for complex shapes. Assume
396                      absolute UTM coordinates.
397                      Format [[x0, y0], [x1, y1], ...]       
398                 
399        Output:       
400            Q: Total flow [m^3/s] across given segments.
401       
402         
403        """       
404       
405        # Find all intersections and associated triangles.
406        segments = self.get_intersecting_segments(polyline, 
407                                                  use_cache=True,
408                                                  verbose=verbose)
409
410        # Get midpoints
411        midpoints = segment_midpoints(segments)       
412       
413        # Make midpoints Geospatial instances
414        midpoints = ensure_geospatial(midpoints, self.geo_reference)       
415       
416        # Compute flow       
417        if verbose: print 'Computing flow through specified cross section'
418       
419        # Get interpolated values
420        xmomentum = self.get_quantity('xmomentum')
421        ymomentum = self.get_quantity('ymomentum')       
422       
423        uh = xmomentum.get_values(interpolation_points=midpoints, use_cache=True)
424        vh = ymomentum.get_values(interpolation_points=midpoints, use_cache=True)
425       
426        # Compute and sum flows across each segment
427        total_flow=0
428        for i in range(len(uh)):
429           
430            # Inner product of momentum vector with segment normal [m^2/s]
431            normal = segments[i].normal
432            normal_momentum = uh[i]*normal[0] + vh[i]*normal[1] 
433               
434            # Flow across this segment [m^3/s]
435            segment_flow = normal_momentum*segments[i].length
436
437            # Accumulate
438            total_flow += segment_flow
439           
440        return total_flow
441
442       
443               
444    def get_energy_through_cross_section(self, polyline,
445                                         kind='total',
446                                         verbose=False):               
447        """Obtain average energy head [m] across specified cross section.
448
449        Inputs:
450            polyline: Representation of desired cross section - it may contain
451                      multiple sections allowing for complex shapes. Assume
452                      absolute UTM coordinates.
453                      Format [[x0, y0], [x1, y1], ...]
454            kind:     Select which energy to compute.
455                      Options are 'specific' and 'total' (default)
456
457        Output:
458            E: Average energy [m] across given segments for all stored times.
459
460        The average velocity is computed for each triangle intersected by
461        the polyline and averaged weighted by segment lengths.
462
463        The typical usage of this function would be to get average energy of
464        flow in a channel, and the polyline would then be a cross section
465        perpendicular to the flow.
466
467        #FIXME (Ole) - need name for this energy reflecting that its dimension
468        is [m].
469        """
470
471        from anuga.config import g, epsilon, velocity_protection as h0       
472                                         
473       
474        # Find all intersections and associated triangles.
475        segments = self.get_intersecting_segments(polyline, 
476                                                  use_cache=True,
477                                                  verbose=verbose)
478
479        # Get midpoints
480        midpoints = segment_midpoints(segments)       
481       
482        # Make midpoints Geospatial instances
483        midpoints = ensure_geospatial(midpoints, self.geo_reference)       
484       
485        # Compute energy       
486        if verbose: print 'Computing %s energy' %kind       
487       
488        # Get interpolated values
489        stage = self.get_quantity('stage')       
490        elevation = self.get_quantity('elevation')               
491        xmomentum = self.get_quantity('xmomentum')
492        ymomentum = self.get_quantity('ymomentum')       
493
494        w = stage.get_values(interpolation_points=midpoints, use_cache=True)
495        z = elevation.get_values(interpolation_points=midpoints, use_cache=True)       
496        uh = xmomentum.get_values(interpolation_points=midpoints, use_cache=True)
497        vh = ymomentum.get_values(interpolation_points=midpoints, use_cache=True)
498        h = w-z # Depth
499       
500        # Compute total length of polyline for use with weighted averages
501        total_line_length = 0.0
502        for segment in segments:
503            total_line_length += segment.length
504       
505        # Compute and sum flows across each segment
506        average_energy=0.0
507        for i in range(len(w)):
508           
509            # Average velocity across this segment
510            if h[i] > epsilon:
511                # Use protection against degenerate velocities
512                u = uh[i]/(h[i] + h0/h[i])
513                v = vh[i]/(h[i] + h0/h[i])
514            else:
515                u = v = 0.0
516               
517            speed_squared = u*u + v*v   
518            kinetic_energy = 0.5*speed_squared/g
519           
520            if kind == 'specific':
521                segment_energy = h[i] + kinetic_energy
522            elif kind == 'total':
523                segment_energy = w[i] + kinetic_energy               
524            else:
525                msg = 'Energy kind must be either "specific" or "total".'
526                msg += ' I got %s' %kind
527               
528
529            # Add to weighted average
530            weigth = segments[i].length/total_line_length
531            average_energy += segment_energy*weigth
532             
533           
534        return average_energy
535       
536
537                       
538
539    def check_integrity(self):
540        Generic_Domain.check_integrity(self)
541
542        #Check that we are solving the shallow water wave equation
543
544        msg = 'First conserved quantity must be "stage"'
545        assert self.conserved_quantities[0] == 'stage', msg
546        msg = 'Second conserved quantity must be "xmomentum"'
547        assert self.conserved_quantities[1] == 'xmomentum', msg
548        msg = 'Third conserved quantity must be "ymomentum"'
549        assert self.conserved_quantities[2] == 'ymomentum', msg
550
551    def extrapolate_second_order_sw(self):
552        #Call correct module function
553        #(either from this module or C-extension)
554        extrapolate_second_order_sw(self)
555
556    def compute_fluxes(self):
557        #Call correct module function
558        #(either from this module or C-extension)
559        compute_fluxes(self)
560
561    def distribute_to_vertices_and_edges(self):
562        # Call correct module function
563        if self.use_edge_limiter:
564            distribute_using_edge_limiter(self)           
565        else:
566            distribute_using_vertex_limiter(self)
567
568
569
570
571    def evolve(self,
572               yieldstep = None,
573               finaltime = None,
574               duration = None,
575               skip_initial_step = False):
576        """Specialisation of basic evolve method from parent class
577        """
578
579        # Call check integrity here rather than from user scripts
580        # self.check_integrity()
581
582        msg = 'Parameter beta_w must be in the interval [0, 2['
583        assert 0 <= self.beta_w <= 2.0, msg
584
585
586        # Initial update of vertex and edge values before any STORAGE
587        # and or visualisation
588        # This is done again in the initialisation of the Generic_Domain
589        # evolve loop but we do it here to ensure the values are ok for storage
590        self.distribute_to_vertices_and_edges()
591
592        if self.store is True and self.time == 0.0:
593            self.initialise_storage()
594            # print 'Storing results in ' + self.writer.filename
595        else:
596            pass
597            # print 'Results will not be stored.'
598            # print 'To store results set domain.store = True'
599            # FIXME: Diagnostic output should be controlled by
600            # a 'verbose' flag living in domain (or in a parent class)
601
602        # Call basic machinery from parent class
603        for t in Generic_Domain.evolve(self,
604                                       yieldstep=yieldstep,
605                                       finaltime=finaltime,
606                                       duration=duration,
607                                       skip_initial_step=skip_initial_step):
608
609            # Store model data, e.g. for subsequent visualisation
610            if self.store is True:
611                self.store_timestep(self.quantities_to_be_stored)
612
613            # FIXME: Could maybe be taken from specified list
614            # of 'store every step' quantities
615
616            # Pass control on to outer loop for more specific actions
617            yield(t)
618
619
620    def initialise_storage(self):
621        """Create and initialise self.writer object for storing data.
622        Also, save x,y and bed elevation
623        """
624
625        from anuga.shallow_water.data_manager import get_dataobject
626
627        # Initialise writer
628        self.writer = get_dataobject(self, mode=netcdf_mode_w)
629
630        # Store vertices and connectivity
631        self.writer.store_connectivity()
632
633
634    def store_timestep(self, name):
635        """Store named quantity and time.
636
637        Precondition:
638           self.write has been initialised
639        """
640        self.writer.store_timestep(name)
641
642       
643    def timestepping_statistics(self,
644                                track_speeds=False,
645                                triangle_id=None):       
646        """Return string with time stepping statistics for printing or logging
647
648        Optional boolean keyword track_speeds decides whether to report
649        location of smallest timestep as well as a histogram and percentile
650        report.
651        """
652
653        from anuga.config import epsilon, g               
654
655
656        # Call basic machinery from parent class
657        msg = Generic_Domain.timestepping_statistics(self,
658                                                     track_speeds,
659                                                     triangle_id)
660
661        if track_speeds is True:
662
663            # qwidth determines the text field used for quantities
664            qwidth = self.qwidth
665       
666            # Selected triangle
667            k = self.k
668
669            # Report some derived quantities at vertices, edges and centroid
670            # specific to the shallow water wave equation
671
672            z = self.quantities['elevation']
673            w = self.quantities['stage']           
674
675            Vw = w.get_values(location='vertices', indices=[k])[0]
676            Ew = w.get_values(location='edges', indices=[k])[0]
677            Cw = w.get_values(location='centroids', indices=[k])
678
679            Vz = z.get_values(location='vertices', indices=[k])[0]
680            Ez = z.get_values(location='edges', indices=[k])[0]
681            Cz = z.get_values(location='centroids', indices=[k])                             
682               
683
684            name = 'depth'
685            Vh = Vw-Vz
686            Eh = Ew-Ez
687            Ch = Cw-Cz
688           
689            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
690                 %(name.ljust(qwidth), Vh[0], Vh[1], Vh[2])
691           
692            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
693                 %(name.ljust(qwidth), Eh[0], Eh[1], Eh[2])
694           
695            s += '    %s: centroid_value = %.4f\n'\
696                 %(name.ljust(qwidth), Ch[0])                               
697           
698            msg += s
699
700            uh = self.quantities['xmomentum']
701            vh = self.quantities['ymomentum']
702
703            Vuh = uh.get_values(location='vertices', indices=[k])[0]
704            Euh = uh.get_values(location='edges', indices=[k])[0]
705            Cuh = uh.get_values(location='centroids', indices=[k])
706           
707            Vvh = vh.get_values(location='vertices', indices=[k])[0]
708            Evh = vh.get_values(location='edges', indices=[k])[0]
709            Cvh = vh.get_values(location='centroids', indices=[k])
710
711            # Speeds in each direction
712            Vu = Vuh/(Vh + epsilon)
713            Eu = Euh/(Eh + epsilon)
714            Cu = Cuh/(Ch + epsilon)           
715            name = 'U'
716            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
717                 %(name.ljust(qwidth), Vu[0], Vu[1], Vu[2])
718           
719            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
720                 %(name.ljust(qwidth), Eu[0], Eu[1], Eu[2])
721           
722            s += '    %s: centroid_value = %.4f\n'\
723                 %(name.ljust(qwidth), Cu[0])                               
724           
725            msg += s
726
727            Vv = Vvh/(Vh + epsilon)
728            Ev = Evh/(Eh + epsilon)
729            Cv = Cvh/(Ch + epsilon)           
730            name = 'V'
731            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
732                 %(name.ljust(qwidth), Vv[0], Vv[1], Vv[2])
733           
734            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
735                 %(name.ljust(qwidth), Ev[0], Ev[1], Ev[2])
736           
737            s += '    %s: centroid_value = %.4f\n'\
738                 %(name.ljust(qwidth), Cv[0])                               
739           
740            msg += s
741
742
743            # Froude number in each direction
744            name = 'Froude (x)'
745            Vfx = Vu/(num.sqrt(g*Vh) + epsilon)
746            Efx = Eu/(num.sqrt(g*Eh) + epsilon)
747            Cfx = Cu/(num.sqrt(g*Ch) + epsilon)
748           
749            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
750                 %(name.ljust(qwidth), Vfx[0], Vfx[1], Vfx[2])
751           
752            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
753                 %(name.ljust(qwidth), Efx[0], Efx[1], Efx[2])
754           
755            s += '    %s: centroid_value = %.4f\n'\
756                 %(name.ljust(qwidth), Cfx[0])                               
757           
758            msg += s
759
760
761            name = 'Froude (y)'
762            Vfy = Vv/(num.sqrt(g*Vh) + epsilon)
763            Efy = Ev/(num.sqrt(g*Eh) + epsilon)
764            Cfy = Cv/(num.sqrt(g*Ch) + epsilon)
765           
766            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
767                 %(name.ljust(qwidth), Vfy[0], Vfy[1], Vfy[2])
768           
769            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
770                 %(name.ljust(qwidth), Efy[0], Efy[1], Efy[2])
771           
772            s += '    %s: centroid_value = %.4f\n'\
773                 %(name.ljust(qwidth), Cfy[0])                               
774           
775            msg += s           
776
777               
778
779        return msg
780       
781       
782
783                def compute_volumetric_balance(self):
784                        """Compute volumetric balance at current timestep.
785                       
786                        Quantities computed are:
787                                Total inflow across boundary
788                                Total outflow across boundary
789                                Total inflow through forcing terms
790                                Total outflow through forcing terms
791                                Current total volume in domain
792                                                               
793                       
794                        """
795               
796                        # Run through boundary array and compute for each segment
797                        # the normal momentum ((uh, vh) dot normal) times segment length.
798                        # Based on sign accumulate this into total_boundary_inflow and total_boundary_outflow.
799                        # The go through explicit forcing update and record the rate of change for stage and
800                        # record into total_forcing_inflow and total_forcing_outflow. Finally compute integral
801                        # of depth to obtain total volume of domain.
802                       
803                       
804                       
805                        pass
806               
807               
808#=============== End of class Shallow Water Domain ===============================
809
810
811#-----------------
812# Flux computation
813#-----------------
814
815def compute_fluxes(domain):
816    """Compute all fluxes and the timestep suitable for all volumes
817    in domain.
818
819    Compute total flux for each conserved quantity using "flux_function"
820
821    Fluxes across each edge are scaled by edgelengths and summed up
822    Resulting flux is then scaled by area and stored in
823    explicit_update for each of the three conserved quantities
824    stage, xmomentum and ymomentum
825
826    The maximal allowable speed computed by the flux_function for each volume
827    is converted to a timestep that must not be exceeded. The minimum of
828    those is computed as the next overall timestep.
829
830    Post conditions:
831      domain.explicit_update is reset to computed flux values
832      domain.timestep is set to the largest step satisfying all volumes.
833   
834
835    This wrapper calls the underlying C version of compute fluxes
836    """
837
838    import sys
839
840    N = len(domain) # number_of_triangles
841
842    # Shortcuts
843    Stage = domain.quantities['stage']
844    Xmom = domain.quantities['xmomentum']
845    Ymom = domain.quantities['ymomentum']
846    Bed = domain.quantities['elevation']
847
848    timestep = float(sys.maxint)
849    from shallow_water_ext import\
850         compute_fluxes_ext_central as compute_fluxes_ext
851
852
853    flux_timestep = compute_fluxes_ext(timestep,
854                                       domain.epsilon,
855                                       domain.H0,
856                                       domain.g,
857                                       domain.neighbours,
858                                       domain.neighbour_edges,
859                                       domain.normals,
860                                       domain.edgelengths,
861                                       domain.radii,
862                                       domain.areas,
863                                       domain.tri_full_flag,
864                                       Stage.edge_values,
865                                       Xmom.edge_values,
866                                       Ymom.edge_values,
867                                       Bed.edge_values,
868                                       Stage.boundary_values,
869                                       Xmom.boundary_values,
870                                       Ymom.boundary_values,
871                                       Stage.explicit_update,
872                                       Xmom.explicit_update,
873                                       Ymom.explicit_update,
874                                       domain.already_computed_flux,
875                                       domain.max_speed,
876                                       int(domain.optimise_dry_cells))
877
878    domain.flux_timestep = flux_timestep
879
880
881
882#---------------------------------------
883# Module functions for gradient limiting
884#---------------------------------------
885
886
887# MH090605 The following method belongs to the shallow_water domain class
888# see comments in the corresponding method in shallow_water_ext.c
889def extrapolate_second_order_sw(domain):
890    """Wrapper calling C version of extrapolate_second_order_sw
891    """
892    import sys
893
894    N = len(domain) # number_of_triangles
895
896    # Shortcuts
897    Stage = domain.quantities['stage']
898    Xmom = domain.quantities['xmomentum']
899    Ymom = domain.quantities['ymomentum']
900    Elevation = domain.quantities['elevation']
901
902    from shallow_water_ext import extrapolate_second_order_sw as extrapol2
903    extrapol2(domain,
904              domain.surrogate_neighbours,
905              domain.number_of_boundaries,
906              domain.centroid_coordinates,
907              Stage.centroid_values,
908              Xmom.centroid_values,
909              Ymom.centroid_values,
910              Elevation.centroid_values,
911              domain.vertex_coordinates,
912              Stage.vertex_values,
913              Xmom.vertex_values,
914              Ymom.vertex_values,
915              Elevation.vertex_values,
916              int(domain.optimise_dry_cells))
917
918
919def distribute_using_vertex_limiter(domain):
920    """Distribution from centroids to vertices specific to the
921    shallow water wave
922    equation.
923
924    It will ensure that h (w-z) is always non-negative even in the
925    presence of steep bed-slopes by taking a weighted average between shallow
926    and deep cases.
927
928    In addition, all conserved quantities get distributed as per either a
929    constant (order==1) or a piecewise linear function (order==2).
930
931    FIXME: more explanation about removal of artificial variability etc
932
933    Precondition:
934      All quantities defined at centroids and bed elevation defined at
935      vertices.
936
937    Postcondition
938      Conserved quantities defined at vertices
939
940    """
941
942   
943
944    # Remove very thin layers of water
945    protect_against_infinitesimal_and_negative_heights(domain)
946
947    # Extrapolate all conserved quantities
948    if domain.optimised_gradient_limiter:
949        # MH090605 if second order,
950        # perform the extrapolation and limiting on
951        # all of the conserved quantities
952
953        if (domain._order_ == 1):
954            for name in domain.conserved_quantities:
955                Q = domain.quantities[name]
956                Q.extrapolate_first_order()
957        elif domain._order_ == 2:
958            domain.extrapolate_second_order_sw()
959        else:
960            raise 'Unknown order'
961    else:
962        # Old code:
963        for name in domain.conserved_quantities:
964            Q = domain.quantities[name]
965
966            if domain._order_ == 1:
967                Q.extrapolate_first_order()
968            elif domain._order_ == 2:
969                Q.extrapolate_second_order_and_limit_by_vertex()
970            else:
971                raise 'Unknown order'
972
973
974    # Take bed elevation into account when water heights are small
975    balance_deep_and_shallow(domain)
976
977    # Compute edge values by interpolation
978    for name in domain.conserved_quantities:
979        Q = domain.quantities[name]
980        Q.interpolate_from_vertices_to_edges()
981
982
983
984def distribute_using_edge_limiter(domain):
985    """Distribution from centroids to edges specific to the
986    shallow water wave
987    equation.
988
989    It will ensure that h (w-z) is always non-negative even in the
990    presence of steep bed-slopes by taking a weighted average between shallow
991    and deep cases.
992
993    In addition, all conserved quantities get distributed as per either a
994    constant (order==1) or a piecewise linear function (order==2).
995
996
997    Precondition:
998      All quantities defined at centroids and bed elevation defined at
999      vertices.
1000
1001    Postcondition
1002      Conserved quantities defined at vertices
1003
1004    """
1005
1006    # Remove very thin layers of water
1007    protect_against_infinitesimal_and_negative_heights(domain)
1008
1009
1010    for name in domain.conserved_quantities:
1011        Q = domain.quantities[name]
1012        if domain._order_ == 1:
1013            Q.extrapolate_first_order()
1014        elif domain._order_ == 2:
1015            Q.extrapolate_second_order_and_limit_by_edge()
1016        else:
1017            raise 'Unknown order'
1018
1019    balance_deep_and_shallow(domain)
1020
1021    # Compute edge values by interpolation
1022    for name in domain.conserved_quantities:
1023        Q = domain.quantities[name]
1024        Q.interpolate_from_vertices_to_edges()
1025
1026
1027def protect_against_infinitesimal_and_negative_heights(domain):
1028    """Protect against infinitesimal heights and associated high velocities
1029    """
1030
1031    # Shortcuts
1032    wc = domain.quantities['stage'].centroid_values
1033    zc = domain.quantities['elevation'].centroid_values
1034    xmomc = domain.quantities['xmomentum'].centroid_values
1035    ymomc = domain.quantities['ymomentum'].centroid_values
1036
1037    from shallow_water_ext import protect
1038
1039    protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
1040            domain.epsilon, wc, zc, xmomc, ymomc)
1041
1042
1043
1044def balance_deep_and_shallow(domain):
1045    """Compute linear combination between stage as computed by
1046    gradient-limiters limiting using w, and stage computed by
1047    gradient-limiters limiting using h (h-limiter).
1048    The former takes precedence when heights are large compared to the
1049    bed slope while the latter takes precedence when heights are
1050    relatively small.  Anything in between is computed as a balanced
1051    linear combination in order to avoid numerical disturbances which
1052    would otherwise appear as a result of hard switching between
1053    modes.
1054
1055    Wrapper for C implementation
1056    """
1057
1058    from shallow_water_ext import balance_deep_and_shallow as balance_deep_and_shallow_c
1059
1060
1061    # Shortcuts
1062    wc = domain.quantities['stage'].centroid_values
1063    zc = domain.quantities['elevation'].centroid_values
1064
1065    wv = domain.quantities['stage'].vertex_values
1066    zv = domain.quantities['elevation'].vertex_values
1067
1068    # Momentums at centroids
1069    xmomc = domain.quantities['xmomentum'].centroid_values
1070    ymomc = domain.quantities['ymomentum'].centroid_values
1071
1072    # Momentums at vertices
1073    xmomv = domain.quantities['xmomentum'].vertex_values
1074    ymomv = domain.quantities['ymomentum'].vertex_values
1075
1076    balance_deep_and_shallow_c(domain,
1077                               wc, zc, wv, zv, wc, 
1078                               xmomc, ymomc, xmomv, ymomv)
1079
1080
1081
1082
1083#------------------------------------------------------------------
1084# Boundary conditions - specific to the shallow water wave equation
1085#------------------------------------------------------------------
1086class Reflective_boundary(Boundary):
1087    """Reflective boundary returns same conserved quantities as
1088    those present in its neighbour volume but reflected.
1089
1090    This class is specific to the shallow water equation as it
1091    works with the momentum quantities assumed to be the second
1092    and third conserved quantities.
1093    """
1094
1095    def __init__(self, domain = None):
1096        Boundary.__init__(self)
1097
1098        if domain is None:
1099            msg = 'Domain must be specified for reflective boundary'
1100            raise msg
1101
1102        # Handy shorthands
1103        self.stage   = domain.quantities['stage'].edge_values
1104        self.xmom    = domain.quantities['xmomentum'].edge_values
1105        self.ymom    = domain.quantities['ymomentum'].edge_values
1106        self.normals = domain.normals
1107
1108        self.conserved_quantities = num.zeros(3, num.Float)
1109
1110    def __repr__(self):
1111        return 'Reflective_boundary'
1112
1113
1114    def evaluate(self, vol_id, edge_id):
1115        """Reflective boundaries reverses the outward momentum
1116        of the volume they serve.
1117        """
1118
1119        q = self.conserved_quantities
1120        q[0] = self.stage[vol_id, edge_id]
1121        q[1] = self.xmom[vol_id, edge_id]
1122        q[2] = self.ymom[vol_id, edge_id]
1123
1124        normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
1125
1126
1127        r = rotate(q, normal, direction = 1)
1128        r[1] = -r[1]
1129        q = rotate(r, normal, direction = -1)
1130
1131        return q
1132
1133
1134
1135
1136class Transmissive_momentum_set_stage_boundary(Boundary):
1137    """Returns same momentum conserved quantities as
1138    those present in its neighbour volume.
1139    Sets stage by specifying a function f of time which may either be a
1140    vector function or a scalar function
1141
1142    Example:
1143
1144    def waveform(t):
1145        return sea_level + normalized_amplitude/cosh(t-25)**2
1146
1147    Bts = Transmissive_momentum_set_stage_boundary(domain, waveform)
1148   
1149
1150    Underlying domain must be specified when boundary is instantiated
1151    """
1152
1153    def __init__(self, domain = None, function=None):
1154        Boundary.__init__(self)
1155
1156        if domain is None:
1157            msg = 'Domain must be specified for this type boundary'
1158            raise msg
1159
1160        if function is None:
1161            msg = 'Function must be specified for this type boundary'
1162            raise msg
1163
1164        self.domain   = domain
1165        self.function = function
1166
1167    def __repr__(self):
1168        return 'Transmissive_momentum_set_stage_boundary(%s)' %self.domain
1169
1170    def evaluate(self, vol_id, edge_id):
1171        """Transmissive momentum set stage boundaries return the edge momentum
1172        values of the volume they serve.
1173        """
1174
1175        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
1176
1177
1178        t = self.domain.get_time()
1179
1180        if hasattr(self.function, 'time'):
1181            # Roll boundary over if time exceeds           
1182            while t > self.function.time[-1]:
1183                msg = 'WARNING: domain time %.2f has exceeded' %t
1184                msg += 'time provided in '
1185                msg += 'transmissive_momentum_set_stage boundary object.\n'
1186                msg += 'I will continue, reusing the object from t==0'
1187                print msg
1188                t -= self.function.time[-1]
1189
1190
1191        value = self.function(t)
1192        try:
1193            x = float(value)
1194        except:
1195            x = float(value[0])
1196           
1197        q[0] = x
1198        return q
1199
1200
1201        # FIXME: Consider this (taken from File_boundary) to allow
1202        # spatial variation
1203        # if vol_id is not None and edge_id is not None:
1204        #     i = self.boundary_indices[ vol_id, edge_id ]
1205        #     return self.F(t, point_id = i)
1206        # else:
1207        #     return self.F(t)
1208
1209
1210# Backward compatibility       
1211# FIXME(Ole): Deprecate
1212class Transmissive_Momentum_Set_Stage_boundary(Transmissive_momentum_set_stage_boundary):
1213    pass
1214
1215     
1216
1217class Transmissive_stage_zero_momentum_boundary(Boundary):
1218    """Return same stage as those present in its neighbour volume. Set momentum to zero.
1219
1220    Underlying domain must be specified when boundary is instantiated
1221    """
1222
1223    def __init__(self, domain=None):
1224        Boundary.__init__(self)
1225
1226        if domain is None:
1227            msg = 'Domain must be specified for '
1228            msg += 'Transmissive_stage_zero_momentum boundary'
1229            raise Exception, msg
1230
1231        self.domain = domain
1232
1233    def __repr__(self):
1234        return 'Transmissive_stage_zero_momentum_boundary(%s)' %self.domain
1235
1236    def evaluate(self, vol_id, edge_id):
1237        """Transmissive boundaries return the edge values
1238        of the volume they serve.
1239        """
1240
1241        q = self.domain.get_conserved_quantities(vol_id, edge=edge_id)
1242       
1243        q[1] = q[2] = 0.0
1244        return q
1245
1246
1247       
1248class Dirichlet_discharge_boundary(Boundary):
1249    """
1250    Sets stage (stage0)
1251    Sets momentum (wh0) in the inward normal direction.
1252
1253    Underlying domain must be specified when boundary is instantiated
1254    """
1255
1256    def __init__(self, domain=None, stage0=None, wh0=None):
1257        Boundary.__init__(self)
1258
1259        if domain is None:
1260            msg = 'Domain must be specified for this type boundary'
1261            raise msg
1262
1263        if stage0 is None:
1264            raise 'set stage'
1265
1266        if wh0 is None:
1267            wh0 = 0.0
1268
1269        self.domain   = domain
1270        self.stage0  = stage0
1271        self.wh0 = wh0
1272
1273    def __repr__(self):
1274        return 'Dirichlet_Discharge_boundary(%s)' %self.domain
1275
1276    def evaluate(self, vol_id, edge_id):
1277        """Set discharge in the (inward) normal direction
1278        """
1279
1280        normal = self.domain.get_normal(vol_id,edge_id)
1281        q = [self.stage0, -self.wh0*normal[0], -self.wh0*normal[1]]
1282        return q
1283
1284
1285        # FIXME: Consider this (taken from File_boundary) to allow
1286        # spatial variation
1287        # if vol_id is not None and edge_id is not None:
1288        #     i = self.boundary_indices[ vol_id, edge_id ]
1289        #     return self.F(t, point_id = i)
1290        # else:
1291        #     return self.F(t)
1292
1293
1294       
1295# Backward compatibility       
1296# FIXME(Ole): Deprecate
1297class Dirichlet_Discharge_boundary(Dirichlet_discharge_boundary):
1298    pass
1299                                                   
1300   
1301
1302       
1303       
1304class Field_boundary(Boundary):
1305    """Set boundary from given field represented in an sww file containing values
1306    for stage, xmomentum and ymomentum.
1307    Optionally, the user can specify mean_stage to offset the stage provided in the
1308    sww file.
1309
1310    This function is a thin wrapper around the generic File_boundary. The
1311    difference between the file_boundary and field_boundary is only that the
1312    field_boundary will allow you to change the level of the stage height when
1313    you read in the boundary condition. This is very useful when running
1314    different tide heights in the same area as you need only to convert one
1315    boundary condition to a SWW file, ideally for tide height of 0 m
1316    (saving disk space). Then you can use field_boundary to read this SWW file
1317    and change the stage height (tide) on the fly depending on the scenario.
1318   
1319    """
1320
1321
1322    def __init__(self, filename, domain,
1323                 mean_stage=0.0,
1324                 time_thinning=1,
1325                 time_limit=None,
1326                 boundary_polygon=None,   
1327                 default_boundary=None,                 
1328                 use_cache=False,
1329                 verbose=False):
1330        """Constructor
1331
1332        filename: Name of sww file
1333        domain: pointer to shallow water domain for which the boundary applies
1334        mean_stage: The mean water level which will be added to stage derived
1335                    from the boundary condition
1336        time_thinning: Will set how many time steps from the sww file read in
1337                       will be interpolated to the boundary. For example if
1338                       the sww file has 1 second time steps and is 24 hours
1339                       in length it has 86400 time steps. If you set
1340                       time_thinning to 1 it will read all these steps.
1341                       If you set it to 100 it will read every 100th step eg
1342                       only 864 step. This parameter is very useful to increase
1343                       the speed of a model run that you are setting up
1344                       and testing.
1345                       
1346        default_boundary: Must be either None or an instance of a
1347                          class descending from class Boundary.
1348                          This will be used in case model time exceeds
1349                          that available in the underlying data.
1350                          Note that mean_stage will also be added to this.
1351                                               
1352        use_cache:
1353        verbose:
1354       
1355        """
1356
1357        # Create generic file_boundary object
1358        self.file_boundary = File_boundary(filename, domain,
1359                                           time_thinning=time_thinning,
1360                                           time_limit=time_limit,
1361                                           boundary_polygon=boundary_polygon,
1362                                           default_boundary=default_boundary,
1363                                           use_cache=use_cache,
1364                                           verbose=verbose)
1365
1366       
1367        # Record information from File_boundary
1368        self.F = self.file_boundary.F
1369        self.domain = self.file_boundary.domain
1370       
1371        # Record mean stage
1372        self.mean_stage = mean_stage
1373
1374
1375    def __repr__(self):
1376        return 'Field boundary'
1377
1378
1379    def evaluate(self, vol_id=None, edge_id=None):
1380        """Return linearly interpolated values based on domain.time
1381
1382        vol_id and edge_id are ignored
1383        """
1384       
1385        # Evaluate file boundary
1386        q = self.file_boundary.evaluate(vol_id, edge_id)
1387
1388        # Adjust stage
1389        for j, name in enumerate(self.domain.conserved_quantities):
1390            if name == 'stage':
1391                q[j] += self.mean_stage
1392        return q
1393
1394   
1395
1396#-----------------------
1397# Standard forcing terms
1398#-----------------------
1399
1400def gravity(domain):
1401    """Apply gravitational pull in the presence of bed slope
1402    Wrapper calls underlying C implementation
1403    """
1404
1405    xmom = domain.quantities['xmomentum'].explicit_update
1406    ymom = domain.quantities['ymomentum'].explicit_update
1407
1408    stage = domain.quantities['stage']
1409    elevation = domain.quantities['elevation']
1410
1411    h = stage.centroid_values - elevation.centroid_values
1412    z = elevation.vertex_values
1413
1414    x = domain.get_vertex_coordinates()
1415    g = domain.g
1416   
1417
1418    from shallow_water_ext import gravity as gravity_c
1419    gravity_c(g, h, z, x, xmom, ymom) #, 1.0e-6)
1420
1421
1422
1423def manning_friction_implicit(domain):
1424    """Apply (Manning) friction to water momentum   
1425    Wrapper for c version
1426    """
1427
1428
1429    #print 'Implicit friction'
1430
1431    xmom = domain.quantities['xmomentum']
1432    ymom = domain.quantities['ymomentum']
1433
1434    w = domain.quantities['stage'].centroid_values
1435    z = domain.quantities['elevation'].centroid_values
1436
1437    uh = xmom.centroid_values
1438    vh = ymom.centroid_values
1439    eta = domain.quantities['friction'].centroid_values
1440
1441    xmom_update = xmom.semi_implicit_update
1442    ymom_update = ymom.semi_implicit_update
1443
1444    N = len(domain)
1445    eps = domain.minimum_allowed_height
1446    g = domain.g
1447
1448    from shallow_water_ext import manning_friction as manning_friction_c
1449    manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
1450
1451
1452def manning_friction_explicit(domain):
1453    """Apply (Manning) friction to water momentum   
1454    Wrapper for c version
1455    """
1456
1457    # print 'Explicit friction'
1458
1459    xmom = domain.quantities['xmomentum']
1460    ymom = domain.quantities['ymomentum']
1461
1462    w = domain.quantities['stage'].centroid_values
1463    z = domain.quantities['elevation'].centroid_values
1464
1465    uh = xmom.centroid_values
1466    vh = ymom.centroid_values
1467    eta = domain.quantities['friction'].centroid_values
1468
1469    xmom_update = xmom.explicit_update
1470    ymom_update = ymom.explicit_update
1471
1472    N = len(domain)
1473    eps = domain.minimum_allowed_height
1474    g = domain.g
1475
1476    from shallow_water_ext import manning_friction as manning_friction_c
1477    manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
1478
1479
1480# FIXME (Ole): This was implemented for use with one of the analytical solutions (Sampson?)
1481# Is it still needed (30 Oct 2007)?
1482def linear_friction(domain):
1483    """Apply linear friction to water momentum
1484
1485    Assumes quantity: 'linear_friction' to be present
1486    """
1487
1488    from math import sqrt
1489
1490    w = domain.quantities['stage'].centroid_values
1491    z = domain.quantities['elevation'].centroid_values
1492    h = w-z
1493
1494    uh = domain.quantities['xmomentum'].centroid_values
1495    vh = domain.quantities['ymomentum'].centroid_values
1496    tau = domain.quantities['linear_friction'].centroid_values
1497
1498    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1499    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1500
1501    N = len(domain) # number_of_triangles
1502    eps = domain.minimum_allowed_height
1503    g = domain.g #Not necessary? Why was this added?
1504
1505    for k in range(N):
1506        if tau[k] >= eps:
1507            if h[k] >= eps:
1508                S = -tau[k]/h[k]
1509
1510                #Update momentum
1511                xmom_update[k] += S*uh[k]
1512                ymom_update[k] += S*vh[k]
1513
1514
1515
1516
1517def depth_dependent_friction(domain, default_friction,
1518                             surface_roughness_data,
1519                             verbose=False):
1520    """Returns an array of friction values for each wet element adjusted for depth.
1521
1522    Inputs:
1523        domain - computational domain object
1524        default_friction - depth independent bottom friction
1525        surface_roughness_data - N x 5 array of n0, d1, n1, d2, n2 values for each
1526        friction region.
1527
1528    Outputs:
1529        wet_friction - Array that can be used directly to update friction as follows:
1530                       domain.set_quantity('friction', wet_friction)
1531
1532       
1533       
1534    """
1535
1536    import Numeric as num
1537   
1538    # Create a temp array to store updated depth dependent friction for wet elements
1539    # EHR this is outwardly inneficient but not obvious how to avoid recreating each call??????
1540    N=len(domain)
1541    wet_friction    = num.zeros(N, num.Float)
1542    wet_friction[:] = default_n0   # Initially assign default_n0 to all array so sure have no zeros values
1543   
1544   
1545    depth = domain.create_quantity_from_expression('stage - elevation')  # create depth instance for this timestep
1546    # Recompute depth as vector 
1547    d = depth.get_values(location='centroids')
1548 
1549    # rebuild the 'friction' values adjusted for depth at this instant
1550    for i in domain.get_wet_elements():                                  # loop for each wet element in domain
1551       
1552        # Get roughness data for each element
1553        n0 = float(surface_roughness_data[i,0])
1554        d1 = float(surface_roughness_data[i,1])
1555        n1 = float(surface_roughness_data[i,2])
1556        d2 = float(surface_roughness_data[i,3])
1557        n2 = float(surface_roughness_data[i,4])
1558       
1559       
1560        # Recompute friction values from depth for this element
1561               
1562        if d[i]   <= d1: 
1563            depth_dependent_friction = n1
1564        elif d[i] >= d2:
1565            depth_dependent_friction = n2
1566        else:
1567            depth_dependent_friction = n1+((n2-n1)/(d2-d1))*(d[i]-d1)
1568           
1569        # check sanity of result
1570        if (depth_dependent_friction  < 0.010 or depth_dependent_friction > 9999.0) :
1571            print model_data.basename+' >>>> WARNING: computed depth_dependent friction out of range ddf,n1,n2 ', depth_dependent_friction, n1,n2
1572       
1573        # update depth dependent friction  for that wet element
1574        wet_friction[i] = depth_dependent_friction
1575       
1576    # EHR add code to show range of 'friction across domain at this instant as sanity check?????????
1577   
1578    if verbose :
1579        nvals=domain.get_quantity('friction').get_values(location='centroids')        # return array of domain nvals
1580        n_min=min(nvals)
1581        n_max=max(nvals)
1582       
1583        print "         ++++ calculate_depth_dependent_friction  - Updated friction - range  %7.3f to %7.3f" %(n_min,n_max)
1584   
1585    return wet_friction
1586
1587
1588
1589
1590
1591
1592#---------------------------------
1593# Experimental auxiliary functions
1594#---------------------------------
1595def check_forcefield(f):
1596    """Check that f is either
1597    1: a callable object f(t,x,y), where x and y are vectors
1598       and that it returns an array or a list of same length
1599       as x and y
1600    2: a scalar
1601    """
1602
1603    if callable(f):
1604        N = 3
1605        x = num.ones(3, num.Float)
1606        y = num.ones(3, num.Float)
1607        try:
1608            q = f(1.0, x=x, y=y)
1609        except Exception, e:
1610            msg = 'Function %s could not be executed:\n%s' %(f, e)
1611            # FIXME: Reconsider this semantics
1612            raise msg
1613
1614        try:
1615            q = num.array(q, num.Float)
1616        except:
1617            msg = 'Return value from vector function %s could ' %f
1618            msg += 'not be converted into a Numeric array of floats.\n'
1619            msg += 'Specified function should return either list or array.'
1620            raise msg
1621
1622        # Is this really what we want?
1623        msg = 'Return vector from function %s ' %f
1624        msg += 'must have same lenght as input vectors'
1625        assert len(q) == N, msg
1626
1627    else:
1628        try:
1629            f = float(f)
1630        except:
1631            msg = 'Force field %s must be either a scalar' %f
1632            msg += ' or a vector function'
1633            raise Exception(msg)
1634    return f
1635
1636
1637class Wind_stress:
1638    """Apply wind stress to water momentum in terms of
1639    wind speed [m/s] and wind direction [degrees]
1640    """
1641
1642    def __init__(self, *args, **kwargs):
1643        """Initialise windfield from wind speed s [m/s]
1644        and wind direction phi [degrees]
1645
1646        Inputs v and phi can be either scalars or Python functions, e.g.
1647
1648        W = Wind_stress(10, 178)
1649
1650        #FIXME - 'normal' degrees are assumed for now, i.e. the
1651        vector (1,0) has zero degrees.
1652        We may need to convert from 'compass' degrees later on and also
1653        map from True north to grid north.
1654
1655        Arguments can also be Python functions of t,x,y as in
1656
1657        def speed(t,x,y):
1658            ...
1659            return s
1660
1661        def angle(t,x,y):
1662            ...
1663            return phi
1664
1665        where x and y are vectors.
1666
1667        and then pass the functions in
1668
1669        W = Wind_stress(speed, angle)
1670
1671        The instantiated object W can be appended to the list of
1672        forcing_terms as in
1673
1674        Alternatively, one vector valued function for (speed, angle)
1675        can be applied, providing both quantities simultaneously.
1676        As in
1677        W = Wind_stress(F), where returns (speed, angle) for each t.
1678
1679        domain.forcing_terms.append(W)
1680        """
1681
1682        from anuga.config import rho_a, rho_w, eta_w
1683
1684        if len(args) == 2:
1685            s = args[0]
1686            phi = args[1]
1687        elif len(args) == 1:
1688            # Assume vector function returning (s, phi)(t,x,y)
1689            vector_function = args[0]
1690            s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
1691            phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
1692        else:
1693           # Assume info is in 2 keyword arguments
1694
1695           if len(kwargs) == 2:
1696               s = kwargs['s']
1697               phi = kwargs['phi']
1698           else:
1699               raise 'Assumes two keyword arguments: s=..., phi=....'
1700
1701        self.speed = check_forcefield(s)
1702        self.phi = check_forcefield(phi)
1703
1704        self.const = eta_w*rho_a/rho_w
1705
1706
1707    def __call__(self, domain):
1708        """Evaluate windfield based on values found in domain
1709        """
1710
1711        from math import pi, cos, sin, sqrt
1712
1713        xmom_update = domain.quantities['xmomentum'].explicit_update
1714        ymom_update = domain.quantities['ymomentum'].explicit_update
1715
1716        N = len(domain) # number_of_triangles
1717        t = domain.time
1718
1719        if callable(self.speed):
1720            xc = domain.get_centroid_coordinates()
1721            s_vec = self.speed(t, xc[:,0], xc[:,1])
1722        else:
1723            # Assume s is a scalar
1724
1725            try:
1726                s_vec = self.speed * num.ones(N, num.Float)
1727            except:
1728                msg = 'Speed must be either callable or a scalar: %s' %self.s
1729                raise msg
1730
1731
1732        if callable(self.phi):
1733            xc = domain.get_centroid_coordinates()
1734            phi_vec = self.phi(t, xc[:,0], xc[:,1])
1735        else:
1736            # Assume phi is a scalar
1737
1738            try:
1739                phi_vec = self.phi * num.ones(N, num.Float)
1740            except:
1741                msg = 'Angle must be either callable or a scalar: %s' %self.phi
1742                raise msg
1743
1744        assign_windfield_values(xmom_update, ymom_update,
1745                                s_vec, phi_vec, self.const)
1746
1747
1748def assign_windfield_values(xmom_update, ymom_update,
1749                            s_vec, phi_vec, const):
1750    """Python version of assigning wind field to update vectors.
1751    A c version also exists (for speed)
1752    """
1753    from math import pi, cos, sin, sqrt
1754
1755    N = len(s_vec)
1756    for k in range(N):
1757        s = s_vec[k]
1758        phi = phi_vec[k]
1759
1760        # Convert to radians
1761        phi = phi*pi/180
1762
1763        # Compute velocity vector (u, v)
1764        u = s*cos(phi)
1765        v = s*sin(phi)
1766
1767        # Compute wind stress
1768        S = const * sqrt(u**2 + v**2)
1769        xmom_update[k] += S*u
1770        ymom_update[k] += S*v
1771
1772
1773
1774
1775
1776class General_forcing:
1777    """General explicit forcing term for update of quantity
1778   
1779    This is used by Inflow and Rainfall for instance
1780   
1781
1782    General_forcing(quantity_name, rate, center, radius, polygon)
1783
1784    domain:     ANUGA computational domain
1785    quantity_name: Name of quantity to update.
1786                   It must be a known conserved quantity.
1787                   
1788    rate [?/s]: Total rate of change over the specified area. 
1789                This parameter can be either a constant or a
1790                function of time. Positive values indicate increases,
1791                negative values indicate decreases.
1792                Rate can be None at initialisation but must be specified
1793                before forcing term is applied (i.e. simulation has started).
1794
1795    center [m]: Coordinates at center of flow point
1796    radius [m]: Size of circular area
1797    polygon:    Arbitrary polygon
1798    default_rate: Rate to be used if rate fails (e.g. if model time exceeds its data)
1799                  Admissible types: None, constant number or function of t
1800
1801
1802    Either center, radius or polygon can be specified but not both.
1803    If neither are specified the entire domain gets updated.
1804    All coordinates to be specified in absolute UTM coordinates (x, y) assuming the zone of domain.   
1805   
1806    Inflow or Rainfall for examples of use
1807    """
1808
1809
1810    # FIXME (AnyOne) : Add various methods to allow spatial variations
1811
1812    def __init__(self,
1813                 domain,
1814                 quantity_name,
1815                 rate=0.0,
1816                 center=None, radius=None,
1817                 polygon=None,
1818                 default_rate=None,
1819                 verbose=False):
1820                     
1821        if center is None:
1822            msg = 'I got radius but no center.'       
1823            assert radius is None, msg
1824           
1825        if radius is None:
1826            msg += 'I got center but no radius.'       
1827            assert center is None, msg
1828           
1829           
1830                     
1831        from math import pi, cos, sin
1832
1833        self.domain = domain
1834        self.quantity_name = quantity_name
1835        self.rate = rate
1836        self.center = ensure_numeric(center)
1837        self.radius = radius
1838        self.polygon = polygon       
1839        self.verbose = verbose
1840        self.value = 0.0 # Can be used to remember value at
1841                         # previous timestep in order to obtain rate
1842
1843        # Get boundary (in absolute coordinates)
1844        bounding_polygon = domain.get_boundary_polygon()
1845
1846
1847        # Update area if applicable
1848        if center is not None and radius is not None:
1849            assert len(center) == 2
1850            msg = 'Polygon cannot be specified when center and radius are'
1851            assert polygon is None, msg
1852
1853
1854            # Check that circle center lies within the mesh.
1855            msg = 'Center %s specified for forcing term did not' %(str(center))
1856            msg += 'fall within the domain boundary.'
1857            assert is_inside_polygon(center, bounding_polygon), msg
1858
1859            # Check that circle periphery lies within the mesh.
1860            N = 100
1861            periphery_points = []
1862            for i in range(N):
1863
1864                theta = 2*pi*i/100
1865               
1866                x = center[0] + radius*cos(theta)
1867                y = center[1] + radius*sin(theta)
1868
1869                periphery_points.append([x,y])
1870
1871
1872            for point in periphery_points:
1873                msg = 'Point %s on periphery for forcing term' %(str(point))
1874                msg += ' did not fall within the domain boundary.'
1875                assert is_inside_polygon(point, bounding_polygon), msg
1876
1877       
1878        if polygon is not None:
1879
1880            # Check that polygon lies within the mesh.
1881            for point in self.polygon:
1882                msg = 'Point %s in polygon for forcing term' %(point)
1883                msg += ' did not fall within the domain boundary.'
1884                assert is_inside_polygon(point, bounding_polygon), msg
1885               
1886                                 
1887
1888        # Pointer to update vector
1889        self.update = domain.quantities[self.quantity_name].explicit_update
1890
1891        # Determine indices in flow area
1892        N = len(domain)   
1893        points = domain.get_centroid_coordinates(absolute=True)
1894
1895        # Calculate indices in exchange area for this forcing term
1896        self.exchange_indices = None
1897        if self.center is not None and self.radius is not None:
1898            # Inlet is circular
1899           
1900            inlet_region = 'center=%s, radius=%s' %(self.center, self.radius) 
1901           
1902            self.exchange_indices = []
1903            for k in range(N):
1904                x, y = points[k,:] # Centroid
1905               
1906                c = self.center
1907                if ((x-c[0])**2+(y-c[1])**2) < self.radius**2:
1908                    self.exchange_indices.append(k)
1909                   
1910        if self.polygon is not None:                   
1911            # Inlet is polygon
1912           
1913            inlet_region = 'polygon=%s' % (self.polygon) 
1914            self.exchange_indices = inside_polygon(points, self.polygon)
1915           
1916           
1917           
1918        if self.exchange_indices is None:
1919            self.exchange_area = polygon_area(bounding_polygon)
1920        else:   
1921            if len(self.exchange_indices) == 0:
1922                msg = 'No triangles have been identified in '
1923                msg += 'specified region: %s' %inlet_region
1924                raise Exception, msg
1925
1926            # Compute exchange area as the sum of areas of triangles identified
1927            # by circle or polygon
1928            self.exchange_area = 0.0
1929            for i in self.exchange_indices:
1930                self.exchange_area += domain.areas[i]
1931           
1932
1933        msg = 'Exchange area in forcing term'
1934        msg += ' has area = %f' %self.exchange_area
1935        assert self.exchange_area > 0.0           
1936           
1937               
1938
1939           
1940        # Check and store default_rate
1941        msg = 'Keyword argument default_rate must be either None ' 
1942        msg += 'or a function of time.\n I got %s' %(str(default_rate)) 
1943        assert default_rate is None or \
1944               type(default_rate) in [IntType, FloatType] or \
1945               callable(default_rate), msg
1946       
1947        if default_rate is not None:
1948
1949            # If it is a constant, make it a function
1950            if not callable(default_rate):
1951                tmp = default_rate
1952                default_rate = lambda t: tmp
1953
1954           
1955            # Check that default_rate is a function of one argument
1956            try:
1957                default_rate(0.0)
1958            except:
1959                raise Exception, msg
1960
1961        self.default_rate = default_rate
1962        self.default_rate_invoked = False    # Flag       
1963       
1964
1965    def __call__(self, domain):
1966        """Apply inflow function at time specified in domain and update stage
1967        """
1968
1969        # Call virtual method allowing local modifications
1970
1971        t = domain.get_time()
1972        try:
1973            rate = self.update_rate(t)
1974        except Modeltime_too_early, e:
1975            raise Modeltime_too_early, e
1976        except Modeltime_too_late, e:
1977            if self.default_rate is None:
1978                raise Exception, e # Reraise exception
1979            else:
1980                # Pass control to default rate function
1981                rate = self.default_rate(t)
1982               
1983                if self.default_rate_invoked is False:
1984                    # Issue warning the first time
1985                    msg = '%s' %str(e)
1986                    msg += 'Instead I will use the default rate: %s\n'\
1987                        %str(self.default_rate) 
1988                    msg += 'Note: Further warnings will be supressed'
1989                    warn(msg)
1990                   
1991                    # FIXME (Ole): Replace this crude flag with
1992                    # Python's ability to print warnings only once.
1993                    # See http://docs.python.org/lib/warning-filter.html
1994                    self.default_rate_invoked = True
1995                   
1996
1997           
1998               
1999
2000        if rate is None:
2001            msg = 'Attribute rate must be specified in General_forcing'
2002            msg += ' or its descendants before attempting to call it'
2003            raise Exception, msg
2004       
2005
2006        # Now rate is a number
2007        if self.verbose is True:
2008            print 'Rate of %s at time = %.2f = %f' %(self.quantity_name,
2009                                                     domain.get_time(),
2010                                                     rate)
2011
2012
2013        if self.exchange_indices is None:
2014            self.update[:] += rate
2015        else:
2016            # Brute force assignment of restricted rate
2017            for k in self.exchange_indices:
2018                self.update[k] += rate
2019
2020
2021    def update_rate(self, t):
2022        """Virtual method allowing local modifications by writing an
2023        overriding version in descendant
2024       
2025        """
2026        if callable(self.rate):
2027            rate = self.rate(t)
2028        else:
2029            rate = self.rate
2030
2031        return rate
2032
2033
2034    def get_quantity_values(self, quantity_name=None):
2035        """Return values for specified quantity restricted to opening
2036       
2037        Optionally a quantity name can be specified if values from another
2038        quantity is sought
2039        """
2040       
2041        if quantity_name is None:
2042            quantity_name = self.quantity_name
2043           
2044        q = self.domain.quantities[quantity_name]
2045        return q.get_values(location='centroids', 
2046                            indices=self.exchange_indices)
2047   
2048
2049    def set_quantity_values(self, val, quantity_name=None):
2050        """Set values for specified quantity restricted to opening
2051       
2052        Optionally a quantity name can be specified if values from another
2053        quantity is sought       
2054        """
2055
2056        if quantity_name is None:
2057            quantity_name = self.quantity_name
2058                   
2059        q = self.domain.quantities[self.quantity_name]               
2060        q.set_values(val, 
2061                     location='centroids', 
2062                     indices=self.exchange_indices)   
2063
2064
2065
2066class Rainfall(General_forcing):
2067    """Class Rainfall - general 'rain over entire domain' forcing term.
2068   
2069    Used for implementing Rainfall over the entire domain.
2070       
2071        Current Limited to only One Gauge..
2072       
2073        Need to add Spatial Varying Capability
2074        (This module came from copying and amending the Inflow Code)
2075   
2076    Rainfall(rain)
2077
2078    domain   
2079    rain [mm/s]:  Total rain rate over the specified domain. 
2080                  NOTE: Raingauge Data needs to reflect the time step.
2081                  IE: if Gauge is mm read at a time step, then the input
2082                  here is as mm/(timeStep) so 10mm in 5minutes becomes
2083                  10/(5x60) = 0.0333mm/s.
2084       
2085       
2086                  This parameter can be either a constant or a
2087                  function of time. Positive values indicate inflow,
2088                  negative values indicate outflow.
2089                  (and be used for Infiltration - Write Seperate Module)
2090                  The specified flow will be divided by the area of
2091                  the inflow region and then applied to update the
2092                  stage quantity.
2093
2094    polygon: Specifies a polygon to restrict the rainfall.
2095   
2096    Examples
2097    How to put them in a run File...
2098       
2099    #------------------------------------------------------------------------
2100    # Setup specialised forcing terms
2101    #------------------------------------------------------------------------
2102    # This is the new element implemented by Ole and Rudy to allow direct
2103    # input of Rainfall in mm/s
2104
2105    catchmentrainfall = Rainfall(rain=file_function('Q100_2hr_Rain.tms')) 
2106                        # Note need path to File in String.
2107                        # Else assumed in same directory
2108
2109    domain.forcing_terms.append(catchmentrainfall)
2110    """
2111
2112   
2113    def __init__(self,
2114                 domain,
2115                 rate=0.0,
2116                 center=None, radius=None,
2117                 polygon=None,
2118                 default_rate=None,                 
2119                 verbose=False):
2120
2121        # Converting mm/s to m/s to apply in ANUGA)
2122        if callable(rate):
2123            rain = lambda t: rate(t)/1000.0
2124        else:
2125            rain = rate/1000.0
2126
2127        if default_rate is not None:   
2128            if callable(default_rate):
2129                default_rain = lambda t: default_rate(t)/1000.0
2130            else:
2131                default_rain = default_rate/1000.0
2132        else:
2133            default_rain = None
2134
2135
2136           
2137           
2138        General_forcing.__init__(self,
2139                                 domain,
2140                                 'stage',
2141                                 rate=rain,
2142                                 center=center, radius=radius,
2143                                 polygon=polygon,
2144                                 default_rate=default_rain,
2145                                 verbose=verbose)
2146
2147       
2148
2149
2150
2151
2152class Inflow(General_forcing):
2153    """Class Inflow - general 'rain and drain' forcing term.
2154   
2155    Useful for implementing flows in and out of the domain.
2156   
2157    Inflow(flow, center, radius, polygon)
2158
2159    domain
2160    rate [m^3/s]: Total flow rate over the specified area. 
2161                  This parameter can be either a constant or a
2162                  function of time. Positive values indicate inflow,
2163                  negative values indicate outflow.
2164                  The specified flow will be divided by the area of
2165                  the inflow region and then applied to update stage.     
2166    center [m]: Coordinates at center of flow point
2167    radius [m]: Size of circular area
2168    polygon:    Arbitrary polygon.
2169
2170    Either center, radius or polygon must be specified
2171   
2172    Examples
2173
2174    # Constant drain at 0.003 m^3/s.
2175    # The outflow area is 0.07**2*pi=0.0154 m^2
2176    # This corresponds to a rate of change of 0.003/0.0154 = 0.2 m/s
2177    #                                     
2178    Inflow((0.7, 0.4), 0.07, -0.003)
2179
2180
2181    # Tap turning up to a maximum inflow of 0.0142 m^3/s.
2182    # The inflow area is 0.03**2*pi = 0.00283 m^2
2183    # This corresponds to a rate of change of 0.0142/0.00283 = 5 m/s     
2184    # over the specified area
2185    Inflow((0.5, 0.5), 0.03, lambda t: min(0.01*t, 0.0142))
2186
2187
2188    #------------------------------------------------------------------------
2189    # Setup specialised forcing terms
2190    #------------------------------------------------------------------------
2191    # This is the new element implemented by Ole to allow direct input
2192    # of Inflow in m^3/s
2193
2194    hydrograph = Inflow(center=(320, 300), radius=10,
2195                        rate=file_function('Q/QPMF_Rot_Sub13.tms'))
2196
2197    domain.forcing_terms.append(hydrograph)
2198   
2199    """
2200
2201
2202    def __init__(self,
2203                 domain,
2204                 rate=0.0,
2205                 center=None, radius=None,
2206                 polygon=None,
2207                 default_rate=None,
2208                 verbose=False):                 
2209
2210
2211        # Create object first to make area is available
2212        General_forcing.__init__(self,
2213                                 domain,
2214                                 'stage',
2215                                 rate=rate,
2216                                 center=center, radius=radius,
2217                                 polygon=polygon,
2218                                 default_rate=default_rate,
2219                                 verbose=verbose)
2220
2221    def update_rate(self, t):
2222        """Virtual method allowing local modifications by writing an
2223        overriding version in descendant
2224
2225        This one converts m^3/s to m/s which can be added directly
2226        to 'stage' in ANUGA
2227        """
2228
2229        if callable(self.rate):
2230            _rate = self.rate(t)/self.exchange_area
2231        else:
2232            _rate = self.rate/self.exchange_area
2233
2234        return _rate
2235
2236
2237
2238
2239#------------------
2240# Initialise module
2241#------------------
2242
2243
2244from anuga.utilities import compile
2245if compile.can_use_C_extension('shallow_water_ext.c'):
2246    # Underlying C implementations can be accessed
2247
2248    from shallow_water_ext import rotate, assign_windfield_values
2249else:
2250    msg = 'C implementations could not be accessed by %s.\n ' %__file__
2251    msg += 'Make sure compile_all.py has been run as described in '
2252    msg += 'the ANUGA installation guide.'
2253    raise Exception, msg
2254
2255
2256# Optimisation with psyco
2257from anuga.config import use_psyco
2258if use_psyco:
2259    try:
2260        import psyco
2261    except:
2262        import os
2263        if os.name == 'posix' and os.uname()[4] in ['x86_64', 'ia64']:
2264            pass
2265            #Psyco isn't supported on 64 bit systems, but it doesn't matter
2266        else:
2267            msg = 'WARNING: psyco (speedup) could not import'+\
2268                  ', you may want to consider installing it'
2269            print msg
2270    else:
2271        psyco.bind(Domain.distribute_to_vertices_and_edges)
2272        psyco.bind(Domain.compute_fluxes)
2273
2274if __name__ == "__main__":
2275    pass
2276
2277
Note: See TracBrowser for help on using the repository browser.