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

Last change on this file since 7352 was 7352, checked in by steve, 15 years ago

Updating cross section

  • Property svn:keywords set to LastChangedDate LastChangedRevision LastChangedBy HeadURL Id
File size: 100.8 KB
Line 
1"""
2Finite-volume computations of the shallow water wave equation.
3
4Title: ANGUA shallow_water_domain - 2D triangular domains for finite-volume
5       computations of the shallow water wave equation.
6
7
8Author: Ole Nielsen, Ole.Nielsen@ga.gov.au
9        Stephen Roberts, Stephen.Roberts@anu.edu.au
10        Duncan Gray, Duncan.Gray@ga.gov.au
11
12CreationDate: 2004
13
14Description:
15    This module contains a specialisation of class Domain from
16    module domain.py consisting of methods specific to the
17    Shallow Water Wave Equation
18
19    U_t + E_x + G_y = S
20
21    where
22
23    U = [w, uh, vh]
24    E = [uh, u^2h + gh^2/2, uvh]
25    G = [vh, uvh, v^2h + gh^2/2]
26    S represents source terms forcing the system
27    (e.g. gravity, friction, wind stress, ...)
28
29    and _t, _x, _y denote the derivative with respect to t, x and y
30    respectively.
31
32
33    The quantities are
34
35    symbol    variable name    explanation
36    x         x                horizontal distance from origin [m]
37    y         y                vertical distance from origin [m]
38    z         elevation        elevation of bed on which flow is modelled [m]
39    h         height           water height above z [m]
40    w         stage            absolute water level, w = z+h [m]
41    u                          speed in the x direction [m/s]
42    v                          speed in the y direction [m/s]
43    uh        xmomentum        momentum in the x direction [m^2/s]
44    vh        ymomentum        momentum in the y direction [m^2/s]
45
46    eta                        mannings friction coefficient [to appear]
47    nu                         wind stress coefficient [to appear]
48
49    The conserved quantities are w, uh, vh
50
51Reference:
52    Catastrophic Collapse of Water Supply Reservoirs in Urban Areas,
53    Christopher Zoppou and Stephen Roberts,
54    Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999
55
56    Hydrodynamic modelling of coastal inundation.
57    Nielsen, O., S. Roberts, D. Gray, A. McPherson and A. Hitchman
58    In Zerger, A. and Argent, R.M. (eds) MODSIM 2005 International Congress on
59    Modelling and Simulation. Modelling and Simulation Society of Australia and
60    New Zealand, December 2005, pp. 518-523. ISBN: 0-9758400-2-9.
61    http://www.mssanz.org.au/modsim05/papers/nielsen.pdf
62
63
64SeeAlso:
65    TRAC administration of ANUGA (User Manuals etc) at
66    https://datamining.anu.edu.au/anuga and Subversion repository at
67    $HeadURL: anuga_core/source/anuga/shallow_water/shallow_water_domain.py $
68
69Constraints: See GPL license in the user guide
70Version: 1.0 ($Revision: 7352 $)
71ModifiedBy:
72    $Author: steve $
73    $Date: 2009-08-11 06:38:46 +0000 (Tue, 11 Aug 2009) $
74"""
75
76# Subversion keywords:
77#
78# $LastChangedDate: 2009-08-11 06:38:46 +0000 (Tue, 11 Aug 2009) $
79# $LastChangedRevision: 7352 $
80# $LastChangedBy: steve $
81
82
83import numpy as num
84
85from anuga.abstract_2d_finite_volumes.neighbour_mesh import segment_midpoints
86from anuga.abstract_2d_finite_volumes.domain import Domain as Generic_Domain
87from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
88     import Boundary
89from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
90     import File_boundary
91from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
92     import Dirichlet_boundary
93from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
94     import Time_boundary
95from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
96     import Transmissive_boundary
97from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
98     import AWI_boundary
99
100from anuga.pmesh.mesh_interface import create_mesh_from_regions
101from anuga.utilities.numerical_tools import gradient, mean, ensure_numeric
102from anuga.geospatial_data.geospatial_data import ensure_geospatial
103
104from anuga.config import minimum_storable_height
105from anuga.config import minimum_allowed_height, maximum_allowed_speed
106from anuga.config import g, epsilon, beta_w, beta_w_dry,\
107     beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, tight_slope_limiters
108from anuga.config import alpha_balance
109from anuga.config import optimise_dry_cells
110from anuga.config import optimised_gradient_limiter
111from anuga.config import use_edge_limiter
112from anuga.config import use_centroid_velocities
113from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
114
115from anuga.fit_interpolate.interpolate import Modeltime_too_late, \
116                                              Modeltime_too_early
117
118from anuga.utilities.polygon import inside_polygon, polygon_area, \
119                                    is_inside_polygon
120import anuga.utilities.log as log
121
122import types
123from types import IntType, FloatType
124from warnings import warn
125
126
127################################################################################
128# Shallow water domain
129################################################################################
130
131##
132# @brief Class for a shallow water domain.
133class Domain(Generic_Domain):
134
135    #conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
136    #other_quantities = ['elevation', 'friction']
137
138    ##
139    # @brief Instantiate a shallow water domain.
140    # @param coordinates
141    # @param vertices
142    # @param boundary
143    # @param tagged_elements
144    # @param geo_reference
145    # @param use_inscribed_circle
146    # @param mesh_filename
147    # @param use_cache
148    # @param verbose
149    # @param full_send_dict
150    # @param ghost_recv_dict
151    # @param processor
152    # @param numproc
153    # @param number_of_full_nodes
154    # @param number_of_full_triangles
155    def __init__(self,
156                 coordinates=None,
157                 vertices=None,
158                 boundary=None,
159                 tagged_elements=None,
160                 geo_reference=None,
161                 use_inscribed_circle=False,
162                 mesh_filename=None,
163                 use_cache=False,
164                 verbose=False,
165                 full_send_dict=None,
166                 ghost_recv_dict=None,
167                 processor=0,
168                 numproc=1,
169                 number_of_full_nodes=None,
170                 number_of_full_triangles=None):
171
172        # Define quantities for the shallow_water domain         
173        conserved_quantities = ['stage', 'xmomentum', 'ymomentum']         
174        other_quantities = ['elevation', 'friction']
175       
176        Generic_Domain.__init__(self,
177                                coordinates,
178                                vertices,
179                                boundary,
180                                conserved_quantities,
181                                other_quantities,
182                                tagged_elements,
183                                geo_reference,
184                                use_inscribed_circle,
185                                mesh_filename,
186                                use_cache,
187                                verbose,
188                                full_send_dict,
189                                ghost_recv_dict,
190                                processor,
191                                numproc,
192                                number_of_full_nodes=number_of_full_nodes,
193                                number_of_full_triangles=number_of_full_triangles)
194
195        self.set_minimum_allowed_height(minimum_allowed_height)
196
197        self.maximum_allowed_speed = maximum_allowed_speed
198        self.g = g
199        self.beta_w = beta_w
200        self.beta_w_dry = beta_w_dry
201        self.beta_uh = beta_uh
202        self.beta_uh_dry = beta_uh_dry
203        self.beta_vh = beta_vh
204        self.beta_vh_dry = beta_vh_dry
205        self.alpha_balance = alpha_balance
206
207        self.tight_slope_limiters = tight_slope_limiters
208        self.optimise_dry_cells = optimise_dry_cells
209
210        self.forcing_terms.append(manning_friction_implicit)
211        self.forcing_terms.append(gravity)
212
213        # Stored output
214        self.store = True
215        self.set_store_vertices_uniquely(False)
216        self.minimum_storable_height = minimum_storable_height
217        self.quantities_to_be_stored = {'elevation': 1, 
218                                        'stage': 2, 
219                                        'xmomentum': 2, 
220                                        'ymomentum': 2}
221
222        # Limiters
223        self.use_edge_limiter = use_edge_limiter
224        self.optimised_gradient_limiter = optimised_gradient_limiter
225        self.use_centroid_velocities = use_centroid_velocities
226
227    ##
228    # @brief
229    # @param beta
230    def set_all_limiters(self, beta):
231        """Shorthand to assign one constant value [0,1] to all limiters.
232        0 Corresponds to first order, where as larger values make use of
233        the second order scheme.
234        """
235
236        self.beta_w = beta
237        self.beta_w_dry = beta
238        self.quantities['stage'].beta = beta
239
240        self.beta_uh = beta
241        self.beta_uh_dry = beta
242        self.quantities['xmomentum'].beta = beta
243
244        self.beta_vh = beta
245        self.beta_vh_dry = beta
246        self.quantities['ymomentum'].beta = beta
247
248    ##
249    # @brief
250    # @param flag
251    # @param reduction
252    def set_store_vertices_uniquely(self, flag, reduction=None):
253        """Decide whether vertex values should be stored uniquely as
254        computed in the model (True) or whether they should be reduced to one
255        value per vertex using self.reduction (False).
256        """
257
258        # FIXME (Ole): how about using the word "continuous vertex values" or
259        # "continuous stage surface"
260        self.smooth = not flag
261
262        # Reduction operation for get_vertex_values
263        if reduction is None:
264            self.reduction = mean
265            #self.reduction = min  #Looks better near steep slopes
266
267    ##
268    # @brief Set the minimum depth that will be written to an SWW file.
269    # @param minimum_storable_height The minimum stored height (in m).
270    def set_minimum_storable_height(self, minimum_storable_height):
271        """Set the minimum depth that will be recognised when writing
272        to an sww file. This is useful for removing thin water layers
273        that seems to be caused by friction creep.
274
275        The minimum allowed sww depth is in meters.
276        """
277
278        self.minimum_storable_height = minimum_storable_height
279
280    ##
281    # @brief
282    # @param minimum_allowed_height
283    def set_minimum_allowed_height(self, minimum_allowed_height):
284        """Set minimum depth that will be recognised in the numerical scheme.
285
286        The minimum allowed depth is in meters.
287
288        The parameter H0 (Minimal height for flux computation)
289        is also set by this function
290        """
291
292        #FIXME (Ole): rename H0 to minimum_allowed_height_in_flux_computation
293
294        #FIXME (Ole): Maybe use histogram to identify isolated extreme speeds
295        #and deal with them adaptively similarly to how we used to use 1 order
296        #steps to recover.
297
298        self.minimum_allowed_height = minimum_allowed_height
299        self.H0 = minimum_allowed_height
300
301    ##
302    # @brief
303    # @param maximum_allowed_speed
304    def set_maximum_allowed_speed(self, maximum_allowed_speed):
305        """Set the maximum particle speed that is allowed in water
306        shallower than minimum_allowed_height. This is useful for
307        controlling speeds in very thin layers of water and at the same time
308        allow some movement avoiding pooling of water.
309        """
310
311        self.maximum_allowed_speed = maximum_allowed_speed
312
313    ##
314    # @brief
315    # @param points_file_block_line_size
316    def set_points_file_block_line_size(self, points_file_block_line_size):
317        """Set the minimum depth that will be recognised when writing
318        to an sww file. This is useful for removing thin water layers
319        that seems to be caused by friction creep.
320
321        The minimum allowed sww depth is in meters.
322        """
323        self.points_file_block_line_size = points_file_block_line_size
324
325       
326    # FIXME: Probably obsolete in its curren form   
327    ##
328    # @brief Set the quantities that will be written to an SWW file.
329    # @param q The quantities to be written.
330    # @note Param 'q' may be None, single quantity or list of quantity strings.
331    # @note If 'q' is None, no quantities will be stored in the SWW file.
332    def set_quantities_to_be_stored(self, q):
333        """Specify which quantities will be stored in the sww file
334       
335        q must be either:
336          - a dictionary with quantity names
337          - a list of quantity names (for backwards compatibility)
338          - None
339
340        The format of the dictionary is as follows
341       
342        quantity_name: flag where flag must be either 1 or 2.
343        If flag is 1, the quantity is considered static and will be
344        stored once at the beginning of the simulation in a 1D array.
345       
346        If flag is 2, the quantity is considered time dependent and
347        it will be stored at each yieldstep by appending it to the
348        appropriate 2D array in the sww file.   
349       
350        If q is None, storage will be switched off altogether.
351       
352        Once the simulation has started and thw sww file opened,
353        this function will have no effect.
354       
355        The format, where q is a list of names is for backwards compatibility only.
356        It will take the specified quantities to be time dependent and assume
357        'elevation' to be static regardless.
358        """
359
360        if q is None:
361            self.quantities_to_be_stored = {}
362            self.store = False
363            return
364
365        # Check correcness
366        for quantity_name in q:
367            msg = ('Quantity %s is not a valid conserved quantity'
368                   % quantity_name)
369            assert quantity_name in self.quantities, msg
370
371        if type(q) == types.ListType:
372
373            msg = 'List arguments to set_quantities_to_be_stored '
374            msg += 'has been deprecated and will be removed in future '
375            msg += 'versions of ANUGA.'
376            msg += 'Please use dictionary argument instead'
377            warn(msg, DeprecationWarning) 
378
379       
380       
381            # FIXME: Raise deprecation warning
382            tmp = {}
383            for x in q:
384                tmp[x] = 2
385            tmp['elevation'] = 1   
386            q = tmp     
387           
388        assert type(q) == types.DictType   
389        self.quantities_to_be_stored = q
390
391    ##
392    # @brief
393    # @param indices
394    def get_wet_elements(self, indices=None):
395        """Return indices for elements where h > minimum_allowed_height
396
397        Optional argument:
398            indices is the set of element ids that the operation applies to.
399
400        Usage:
401            indices = get_wet_elements()
402
403        Note, centroid values are used for this operation
404        """
405
406        # Water depth below which it is considered to be 0 in the model
407        # FIXME (Ole): Allow this to be specified as a keyword argument as well
408        from anuga.config import minimum_allowed_height
409
410        elevation = self.get_quantity('elevation').\
411                        get_values(location='centroids', indices=indices)
412        stage = self.get_quantity('stage').\
413                    get_values(location='centroids', indices=indices)
414        depth = stage - elevation
415
416        # Select indices for which depth > 0
417        wet_indices = num.compress(depth > minimum_allowed_height,
418                                   num.arange(len(depth)))
419        return wet_indices
420
421    ##
422    # @brief
423    # @param indices
424    def get_maximum_inundation_elevation(self, indices=None):
425        """Return highest elevation where h > 0
426
427        Optional argument:
428            indices is the set of element ids that the operation applies to.
429
430        Usage:
431            q = get_maximum_inundation_elevation()
432
433        Note, centroid values are used for this operation
434        """
435
436        wet_elements = self.get_wet_elements(indices)
437        return self.get_quantity('elevation').\
438                   get_maximum_value(indices=wet_elements)
439
440    ##
441    # @brief
442    # @param indices
443    def get_maximum_inundation_location(self, indices=None):
444        """Return location of highest elevation where h > 0
445
446        Optional argument:
447            indices is the set of element ids that the operation applies to.
448
449        Usage:
450            q = get_maximum_inundation_location()
451
452        Note, centroid values are used for this operation
453        """
454
455        wet_elements = self.get_wet_elements(indices)
456        return self.get_quantity('elevation').\
457                   get_maximum_location(indices=wet_elements)
458
459
460
461    ##
462    # @brief Get the total flow through an arbitrary poly line.
463    # @param polyline Representation of desired cross section.
464    # @param verbose True if this method is to be verbose.
465    # @note 'polyline' may contain multiple sections allowing complex shapes.
466    # @note Assume absolute UTM coordinates.
467    def get_flow_through_cross_section(self, polyline, verbose=False):
468        """Get the total flow through an arbitrary poly line.
469
470        This is a run-time equivalent of the function with same name
471        in data_manager.py
472
473        Input:
474            polyline: Representation of desired cross section - it may contain
475                      multiple sections allowing for complex shapes. Assume
476                      absolute UTM coordinates.
477                      Format [[x0, y0], [x1, y1], ...]
478
479        Output:
480            Q: Total flow [m^3/s] across given segments.
481        """
482
483
484        cross_section = Cross_section(self, polyline, verbose)
485
486        return cross_section.get_flow_through_cross_section()
487
488
489
490
491
492    ##
493    # @brief Get the total flow through an arbitrary poly line.
494    # @param polyline Representation of desired cross section.
495    # @param verbose True if this method is to be verbose.
496    # @note 'polyline' may contain multiple sections allowing complex shapes.
497    # @note Assume absolute UTM coordinates.
498    def old_get_flow_through_cross_section(self, polyline, verbose=False):
499        """Get the total flow through an arbitrary poly line.
500
501        This is a run-time equivalent of the function with same name
502        in data_manager.py
503
504        Input:
505            polyline: Representation of desired cross section - it may contain
506                      multiple sections allowing for complex shapes. Assume
507                      absolute UTM coordinates.
508                      Format [[x0, y0], [x1, y1], ...]
509
510        Output:
511            Q: Total flow [m^3/s] across given segments.
512        """
513
514        # Find all intersections and associated triangles.
515        segments = self.get_intersecting_segments(polyline, use_cache=True,
516                                                  verbose=verbose)
517
518        # Get midpoints
519        midpoints = segment_midpoints(segments)
520
521        # Make midpoints Geospatial instances
522        midpoints = ensure_geospatial(midpoints, self.geo_reference)
523
524        # Compute flow
525        if verbose:
526            log.critical('Computing flow through specified cross section')
527
528        # Get interpolated values
529        xmomentum = self.get_quantity('xmomentum')
530        ymomentum = self.get_quantity('ymomentum')
531
532        uh = xmomentum.get_values(interpolation_points=midpoints,
533                                  use_cache=True)
534        vh = ymomentum.get_values(interpolation_points=midpoints,
535                                  use_cache=True)
536
537        # Compute and sum flows across each segment
538        total_flow = 0
539        for i in range(len(uh)):
540            # Inner product of momentum vector with segment normal [m^2/s]
541            normal = segments[i].normal
542            normal_momentum = uh[i]*normal[0] + vh[i]*normal[1]
543
544            # Flow across this segment [m^3/s]
545            segment_flow = normal_momentum*segments[i].length
546
547            # Accumulate
548            total_flow += segment_flow
549
550        return total_flow
551
552    ##
553    # @brief
554    # @param polyline Representation of desired cross section.
555    # @param kind Select energy type to compute ('specific' or 'total').
556    # @param verbose True if this method is to be verbose.
557    # @note 'polyline' may contain multiple sections allowing complex shapes.
558    # @note Assume absolute UTM coordinates.
559    def new_get_energy_through_cross_section(self, polyline,
560                                         kind='total',
561                                         verbose=False):
562        """Obtain average energy head [m] across specified cross section.
563
564        Inputs:
565            polyline: Representation of desired cross section - it may contain
566                      multiple sections allowing for complex shapes. Assume
567                      absolute UTM coordinates.
568                      Format [[x0, y0], [x1, y1], ...]
569            kind:     Select which energy to compute.
570                      Options are 'specific' and 'total' (default)
571
572        Output:
573            E: Average energy [m] across given segments for all stored times.
574
575        The average velocity is computed for each triangle intersected by
576        the polyline and averaged weighted by segment lengths.
577
578        The typical usage of this function would be to get average energy of
579        flow in a channel, and the polyline would then be a cross section
580        perpendicular to the flow.
581
582        #FIXME (Ole) - need name for this energy reflecting that its dimension
583        is [m].
584        """
585
586
587
588        cross_section = Cross_section(self, polyline, verbose)
589
590        return cross_section.get_energy_through_cross_section()
591
592
593    ##
594    # @brief
595    # @param polyline Representation of desired cross section.
596    # @param kind Select energy type to compute ('specific' or 'total').
597    # @param verbose True if this method is to be verbose.
598    # @note 'polyline' may contain multiple sections allowing complex shapes.
599    # @note Assume absolute UTM coordinates.
600    def get_energy_through_cross_section(self, polyline,
601                                         kind='total',
602                                         verbose=False):
603        """Obtain average energy head [m] across specified cross section.
604
605        Inputs:
606            polyline: Representation of desired cross section - it may contain
607                      multiple sections allowing for complex shapes. Assume
608                      absolute UTM coordinates.
609                      Format [[x0, y0], [x1, y1], ...]
610            kind:     Select which energy to compute.
611                      Options are 'specific' and 'total' (default)
612
613        Output:
614            E: Average energy [m] across given segments for all stored times.
615
616        The average velocity is computed for each triangle intersected by
617        the polyline and averaged weighted by segment lengths.
618
619        The typical usage of this function would be to get average energy of
620        flow in a channel, and the polyline would then be a cross section
621        perpendicular to the flow.
622
623        #FIXME (Ole) - need name for this energy reflecting that its dimension
624        is [m].
625        """
626
627        from anuga.config import g, epsilon, velocity_protection as h0
628
629        # Find all intersections and associated triangles.
630        segments = self.get_intersecting_segments(polyline, use_cache=True,
631                                                  verbose=verbose)
632
633        # Get midpoints
634        midpoints = segment_midpoints(segments)
635
636        # Make midpoints Geospatial instances
637        midpoints = ensure_geospatial(midpoints, self.geo_reference)
638
639        # Compute energy
640        if verbose: log.critical('Computing %s energy' % kind)
641
642        # Get interpolated values
643        stage = self.get_quantity('stage')
644        elevation = self.get_quantity('elevation')
645        xmomentum = self.get_quantity('xmomentum')
646        ymomentum = self.get_quantity('ymomentum')
647
648        w = stage.get_values(interpolation_points=midpoints, use_cache=True)
649        z = elevation.get_values(interpolation_points=midpoints, use_cache=True)
650        uh = xmomentum.get_values(interpolation_points=midpoints,
651                                  use_cache=True)
652        vh = ymomentum.get_values(interpolation_points=midpoints,
653                                  use_cache=True)
654        h = w-z                # Depth
655
656        # Compute total length of polyline for use with weighted averages
657        total_line_length = 0.0
658        for segment in segments:
659            total_line_length += segment.length
660
661        # Compute and sum flows across each segment
662        average_energy = 0.0
663        for i in range(len(w)):
664            # Average velocity across this segment
665            if h[i] > epsilon:
666                # Use protection against degenerate velocities
667                u = uh[i]/(h[i] + h0/h[i])
668                v = vh[i]/(h[i] + h0/h[i])
669            else:
670                u = v = 0.0
671
672            speed_squared = u*u + v*v
673            kinetic_energy = 0.5*speed_squared/g
674
675            if kind == 'specific':
676                segment_energy = h[i] + kinetic_energy
677            elif kind == 'total':
678                segment_energy = w[i] + kinetic_energy
679            else:
680                msg = 'Energy kind must be either "specific" or "total".'
681                msg += ' I got %s' %kind
682
683            # Add to weighted average
684            weigth = segments[i].length/total_line_length
685            average_energy += segment_energy*weigth
686
687        return average_energy
688
689
690    ##
691    # @brief Run integrity checks on shallow water domain.
692    def check_integrity(self):
693        Generic_Domain.check_integrity(self)
694
695        #Check that we are solving the shallow water wave equation
696        msg = 'First conserved quantity must be "stage"'
697        assert self.conserved_quantities[0] == 'stage', msg
698        msg = 'Second conserved quantity must be "xmomentum"'
699        assert self.conserved_quantities[1] == 'xmomentum', msg
700        msg = 'Third conserved quantity must be "ymomentum"'
701        assert self.conserved_quantities[2] == 'ymomentum', msg
702
703    ##
704    # @brief
705    def extrapolate_second_order_sw(self):
706        #Call correct module function (either from this module or C-extension)
707        extrapolate_second_order_sw(self)
708
709    ##
710    # @brief
711    def compute_fluxes(self):
712        #Call correct module function (either from this module or C-extension)
713        compute_fluxes(self)
714
715    ##
716    # @brief
717    def distribute_to_vertices_and_edges(self):
718        # Call correct module function
719        if self.use_edge_limiter:
720            distribute_using_edge_limiter(self)
721        else:
722            distribute_using_vertex_limiter(self)
723
724
725
726    ##
727    # @brief Evolve the model by one step.
728    # @param yieldstep
729    # @param finaltime
730    # @param duration
731    # @param skip_initial_step
732    def evolve(self,
733               yieldstep=None,
734               finaltime=None,
735               duration=None,
736               skip_initial_step=False):
737        """Specialisation of basic evolve method from parent class"""
738
739        # Call check integrity here rather than from user scripts
740        # self.check_integrity()
741
742        msg = 'Attribute self.beta_w must be in the interval [0, 2]'
743        assert 0 <= self.beta_w <= 2.0, msg
744
745        # Initial update of vertex and edge values before any STORAGE
746        # and or visualisation.
747        # This is done again in the initialisation of the Generic_Domain
748        # evolve loop but we do it here to ensure the values are ok for storage.
749        self.distribute_to_vertices_and_edges()
750
751        if self.store is True and self.time == 0.0:
752            self.initialise_storage()
753
754        # Call basic machinery from parent class
755        for t in Generic_Domain.evolve(self, yieldstep=yieldstep,
756                                       finaltime=finaltime, duration=duration,
757                                       skip_initial_step=skip_initial_step):
758            # Store model data, e.g. for subsequent visualisation
759            if self.store is True:
760                self.store_timestep()
761
762            # Pass control on to outer loop for more specific actions
763            yield(t)
764
765    ##
766    # @brief
767    def initialise_storage(self):
768        """Create and initialise self.writer object for storing data.
769        Also, save x,y and bed elevation
770        """
771
772        from anuga.shallow_water.data_manager import SWW_file
773       
774        # Initialise writer
775        self.writer = SWW_file(self)
776
777        # Store vertices and connectivity
778        self.writer.store_connectivity()
779
780    ##
781    # @brief
782    # @param name
783    def store_timestep(self):
784        """Store time dependent quantities and time.
785
786        Precondition:
787           self.writer has been initialised
788        """
789
790        self.writer.store_timestep()
791
792    ##
793    # @brief Get time stepping statistics string for printing.
794    # @param track_speeds
795    # @param triangle_id
796    def timestepping_statistics(self,
797                                track_speeds=False,
798                                triangle_id=None):
799        """Return string with time stepping statistics for printing or logging
800
801        Optional boolean keyword track_speeds decides whether to report
802        location of smallest timestep as well as a histogram and percentile
803        report.
804        """
805
806        from anuga.config import epsilon, g
807
808        # Call basic machinery from parent class
809        msg = Generic_Domain.timestepping_statistics(self, track_speeds,
810                                                     triangle_id)
811
812        if track_speeds is True:
813            # qwidth determines the text field used for quantities
814            qwidth = self.qwidth
815
816            # Selected triangle
817            k = self.k
818
819            # Report some derived quantities at vertices, edges and centroid
820            # specific to the shallow water wave equation
821            z = self.quantities['elevation']
822            w = self.quantities['stage']
823
824            Vw = w.get_values(location='vertices', indices=[k])[0]
825            Ew = w.get_values(location='edges', indices=[k])[0]
826            Cw = w.get_values(location='centroids', indices=[k])
827
828            Vz = z.get_values(location='vertices', indices=[k])[0]
829            Ez = z.get_values(location='edges', indices=[k])[0]
830            Cz = z.get_values(location='centroids', indices=[k])
831
832            name = 'depth'
833            Vh = Vw-Vz
834            Eh = Ew-Ez
835            Ch = Cw-Cz
836
837            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
838                 %(name.ljust(qwidth), Vh[0], Vh[1], Vh[2])
839
840            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
841                 %(name.ljust(qwidth), Eh[0], Eh[1], Eh[2])
842
843            s += '    %s: centroid_value = %.4f\n'\
844                 %(name.ljust(qwidth), Ch[0])
845
846            msg += s
847
848            uh = self.quantities['xmomentum']
849            vh = self.quantities['ymomentum']
850
851            Vuh = uh.get_values(location='vertices', indices=[k])[0]
852            Euh = uh.get_values(location='edges', indices=[k])[0]
853            Cuh = uh.get_values(location='centroids', indices=[k])
854
855            Vvh = vh.get_values(location='vertices', indices=[k])[0]
856            Evh = vh.get_values(location='edges', indices=[k])[0]
857            Cvh = vh.get_values(location='centroids', indices=[k])
858
859            # Speeds in each direction
860            Vu = Vuh/(Vh + epsilon)
861            Eu = Euh/(Eh + epsilon)
862            Cu = Cuh/(Ch + epsilon)
863            name = 'U'
864            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
865                 %(name.ljust(qwidth), Vu[0], Vu[1], Vu[2])
866
867            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
868                 %(name.ljust(qwidth), Eu[0], Eu[1], Eu[2])
869
870            s += '    %s: centroid_value = %.4f\n'\
871                 %(name.ljust(qwidth), Cu[0])
872
873            msg += s
874
875            Vv = Vvh/(Vh + epsilon)
876            Ev = Evh/(Eh + epsilon)
877            Cv = Cvh/(Ch + epsilon)
878            name = 'V'
879            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
880                 %(name.ljust(qwidth), Vv[0], Vv[1], Vv[2])
881
882            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
883                 %(name.ljust(qwidth), Ev[0], Ev[1], Ev[2])
884
885            s += '    %s: centroid_value = %.4f\n'\
886                 %(name.ljust(qwidth), Cv[0])
887
888            msg += s
889
890            # Froude number in each direction
891            name = 'Froude (x)'
892            Vfx = Vu/(num.sqrt(g*Vh) + epsilon)
893            Efx = Eu/(num.sqrt(g*Eh) + epsilon)
894            Cfx = Cu/(num.sqrt(g*Ch) + epsilon)
895
896            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
897                 %(name.ljust(qwidth), Vfx[0], Vfx[1], Vfx[2])
898
899            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
900                 %(name.ljust(qwidth), Efx[0], Efx[1], Efx[2])
901
902            s += '    %s: centroid_value = %.4f\n'\
903                 %(name.ljust(qwidth), Cfx[0])
904
905            msg += s
906
907            name = 'Froude (y)'
908            Vfy = Vv/(num.sqrt(g*Vh) + epsilon)
909            Efy = Ev/(num.sqrt(g*Eh) + epsilon)
910            Cfy = Cv/(num.sqrt(g*Ch) + epsilon)
911
912            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
913                 %(name.ljust(qwidth), Vfy[0], Vfy[1], Vfy[2])
914
915            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
916                 %(name.ljust(qwidth), Efy[0], Efy[1], Efy[2])
917
918            s += '    %s: centroid_value = %.4f\n'\
919                 %(name.ljust(qwidth), Cfy[0])
920
921            msg += s
922
923        return msg
924       
925       
926
927    def compute_boundary_flows(self):
928        """Compute boundary flows at current timestep.
929                       
930        Quantities computed are:
931           Total inflow across boundary
932           Total outflow across boundary
933           Flow across each tagged boundary segment
934        """
935               
936        # Run through boundary array and compute for each segment
937        # the normal momentum ((uh, vh) dot normal) times segment length.
938        # Based on sign accumulate this into boundary_inflow and boundary_outflow.
939                       
940        # Compute flows along boundary
941       
942        uh = self.get_quantity('xmomentum').get_values(location='edges')
943        vh = self.get_quantity('ymomentum').get_values(location='edges')       
944       
945        # Loop through edges that lie on the boundary and calculate
946        # flows
947        boundary_flows = {}
948        total_boundary_inflow = 0.0
949        total_boundary_outflow = 0.0
950        for vol_id, edge_id in self.boundary:
951            # Compute normal flow across edge. Since normal vector points
952            # away from triangle, a positive sign means that water
953            # flows *out* from this triangle.
954             
955            momentum = [uh[vol_id, edge_id], vh[vol_id, edge_id]]
956            normal = self.mesh.get_normal(vol_id, edge_id)
957            length = self.mesh.get_edgelength(vol_id, edge_id)           
958            normal_flow = num.dot(momentum, normal)*length
959           
960            # Reverse sign so that + is taken to mean inflow
961            # and - means outflow. This is more intuitive.
962            edge_flow = -normal_flow
963           
964            # Tally up inflows and outflows separately
965            if edge_flow > 0:
966                # Flow is inflow     
967                total_boundary_inflow += edge_flow                                 
968            else:
969                # Flow is outflow
970                total_boundary_outflow += edge_flow   
971
972            # Tally up flows by boundary tag
973            tag = self.boundary[(vol_id, edge_id)]
974           
975            if tag not in boundary_flows:
976                boundary_flows[tag] = 0.0
977            boundary_flows[tag] += edge_flow
978           
979               
980        return boundary_flows, total_boundary_inflow, total_boundary_outflow
981       
982
983    def compute_forcing_flows(self):
984        """Compute flows in and out of domain due to forcing terms.
985                       
986        Quantities computed are:
987               
988       
989           Total inflow through forcing terms
990           Total outflow through forcing terms
991           Current total volume in domain       
992
993        """
994
995        #FIXME(Ole): We need to separate what part of explicit_update was
996        # due to the normal flux calculations and what is due to forcing terms.
997       
998        pass
999                       
1000       
1001    def compute_total_volume(self):
1002        """Compute total volume (m^3) of water in entire domain
1003        """
1004       
1005        area = self.mesh.get_areas()
1006        volume = 0.0
1007       
1008        stage = self.get_quantity('stage').get_values(location='centroids')
1009        elevation = self.get_quantity('elevation').get_values(location='centroids')       
1010        depth = stage-elevation
1011       
1012        return num.sum(depth*area)
1013       
1014       
1015    def volumetric_balance_statistics(self):               
1016        """Create volumetric balance report suitable for printing or logging.
1017        """
1018       
1019        (boundary_flows, total_boundary_inflow,
1020         total_boundary_outflow) = self.compute_boundary_flows() 
1021       
1022        s = '---------------------------\n'       
1023        s += 'Volumetric balance report:\n'
1024        s += '--------------------------\n'
1025        s += 'Total boundary inflow [m^3/s]: %.2f\n' % total_boundary_inflow
1026        s += 'Total boundary outflow [m^3/s]: %.2f\n' % total_boundary_outflow       
1027        s += 'Net boundary flow by tags [m^3/s]\n'
1028        for tag in boundary_flows:
1029            s += '    %s [m^3/s]: %.2f\n' % (tag, boundary_flows[tag])
1030       
1031        s += 'Total net boundary flow [m^3/s]: %.2f\n' % (total_boundary_inflow + total_boundary_outflow) 
1032        s += 'Total volume in domain [m^3]: %.2f\n' % self.compute_total_volume()
1033       
1034        # The go through explicit forcing update and record the rate of change for stage and
1035        # record into forcing_inflow and forcing_outflow. Finally compute integral
1036        # of depth to obtain total volume of domain.
1037       
1038        # FIXME(Ole): This part is not yet done.               
1039       
1040        return s       
1041           
1042################################################################################
1043# End of class Shallow Water Domain
1044################################################################################
1045
1046#-----------------
1047# Flux computation
1048#-----------------
1049
1050## @brief Compute fluxes and timestep suitable for all volumes in domain.
1051# @param domain The domain to calculate fluxes for.
1052def compute_fluxes(domain):
1053    """Compute fluxes and timestep suitable for all volumes in domain.
1054
1055    Compute total flux for each conserved quantity using "flux_function"
1056
1057    Fluxes across each edge are scaled by edgelengths and summed up
1058    Resulting flux is then scaled by area and stored in
1059    explicit_update for each of the three conserved quantities
1060    stage, xmomentum and ymomentum
1061
1062    The maximal allowable speed computed by the flux_function for each volume
1063    is converted to a timestep that must not be exceeded. The minimum of
1064    those is computed as the next overall timestep.
1065
1066    Post conditions:
1067      domain.explicit_update is reset to computed flux values
1068      domain.timestep is set to the largest step satisfying all volumes.
1069
1070    This wrapper calls the underlying C version of compute fluxes
1071    """
1072
1073    import sys
1074    from shallow_water_ext import compute_fluxes_ext_central \
1075                                  as compute_fluxes_ext
1076
1077    N = len(domain)    # number_of_triangles
1078
1079    # Shortcuts
1080    Stage = domain.quantities['stage']
1081    Xmom = domain.quantities['xmomentum']
1082    Ymom = domain.quantities['ymomentum']
1083    Bed = domain.quantities['elevation']
1084
1085    timestep = float(sys.maxint)
1086
1087    flux_timestep = compute_fluxes_ext(timestep,
1088                                       domain.epsilon,
1089                                       domain.H0,
1090                                       domain.g,
1091                                       domain.neighbours,
1092                                       domain.neighbour_edges,
1093                                       domain.normals,
1094                                       domain.edgelengths,
1095                                       domain.radii,
1096                                       domain.areas,
1097                                       domain.tri_full_flag,
1098                                       Stage.edge_values,
1099                                       Xmom.edge_values,
1100                                       Ymom.edge_values,
1101                                       Bed.edge_values,
1102                                       Stage.boundary_values,
1103                                       Xmom.boundary_values,
1104                                       Ymom.boundary_values,
1105                                       Stage.explicit_update,
1106                                       Xmom.explicit_update,
1107                                       Ymom.explicit_update,
1108                                       domain.already_computed_flux,
1109                                       domain.max_speed,
1110                                       int(domain.optimise_dry_cells))
1111
1112    domain.flux_timestep = flux_timestep
1113
1114################################################################################
1115# Module functions for gradient limiting
1116################################################################################
1117
1118##
1119# @brief Wrapper for C version of extrapolate_second_order_sw.
1120# @param domain The domain to operate on.
1121# @note MH090605 The following method belongs to the shallow_water domain class
1122#       see comments in the corresponding method in shallow_water_ext.c
1123def extrapolate_second_order_sw(domain):
1124    """Wrapper calling C version of extrapolate_second_order_sw"""
1125
1126    import sys
1127    from shallow_water_ext import extrapolate_second_order_sw as extrapol2
1128
1129    N = len(domain) # number_of_triangles
1130
1131    # Shortcuts
1132    Stage = domain.quantities['stage']
1133    Xmom = domain.quantities['xmomentum']
1134    Ymom = domain.quantities['ymomentum']
1135    Elevation = domain.quantities['elevation']
1136
1137    extrapol2(domain,
1138              domain.surrogate_neighbours,
1139              domain.number_of_boundaries,
1140              domain.centroid_coordinates,
1141              Stage.centroid_values,
1142              Xmom.centroid_values,
1143              Ymom.centroid_values,
1144              Elevation.centroid_values,
1145              domain.vertex_coordinates,
1146              Stage.vertex_values,
1147              Xmom.vertex_values,
1148              Ymom.vertex_values,
1149              Elevation.vertex_values,
1150              int(domain.optimise_dry_cells))
1151
1152##
1153# @brief Distribution from centroids to vertices specific to the SWW eqn.
1154# @param domain The domain to operate on.
1155def distribute_using_vertex_limiter(domain):
1156    """Distribution from centroids to vertices specific to the SWW equation.
1157
1158    It will ensure that h (w-z) is always non-negative even in the
1159    presence of steep bed-slopes by taking a weighted average between shallow
1160    and deep cases.
1161
1162    In addition, all conserved quantities get distributed as per either a
1163    constant (order==1) or a piecewise linear function (order==2).
1164
1165    FIXME: more explanation about removal of artificial variability etc
1166
1167    Precondition:
1168      All quantities defined at centroids and bed elevation defined at
1169      vertices.
1170
1171    Postcondition
1172      Conserved quantities defined at vertices
1173    """
1174
1175    # Remove very thin layers of water
1176    protect_against_infinitesimal_and_negative_heights(domain)
1177
1178    # Extrapolate all conserved quantities
1179    if domain.optimised_gradient_limiter:
1180        # MH090605 if second order,
1181        # perform the extrapolation and limiting on
1182        # all of the conserved quantities
1183
1184        if (domain._order_ == 1):
1185            for name in domain.conserved_quantities:
1186                Q = domain.quantities[name]
1187                Q.extrapolate_first_order()
1188        elif domain._order_ == 2:
1189            domain.extrapolate_second_order_sw()
1190        else:
1191            raise 'Unknown order'
1192    else:
1193        # Old code:
1194        for name in domain.conserved_quantities:
1195            Q = domain.quantities[name]
1196
1197            if domain._order_ == 1:
1198                Q.extrapolate_first_order()
1199            elif domain._order_ == 2:
1200                Q.extrapolate_second_order_and_limit_by_vertex()
1201            else:
1202                raise 'Unknown order'
1203
1204    # Take bed elevation into account when water heights are small
1205    balance_deep_and_shallow(domain)
1206
1207    # Compute edge values by interpolation
1208    for name in domain.conserved_quantities:
1209        Q = domain.quantities[name]
1210        Q.interpolate_from_vertices_to_edges()
1211
1212##
1213# @brief Distribution from centroids to edges specific to the SWW eqn.
1214# @param domain The domain to operate on.
1215def distribute_using_edge_limiter(domain):
1216    """Distribution from centroids to edges specific to the SWW eqn.
1217
1218    It will ensure that h (w-z) is always non-negative even in the
1219    presence of steep bed-slopes by taking a weighted average between shallow
1220    and deep cases.
1221
1222    In addition, all conserved quantities get distributed as per either a
1223    constant (order==1) or a piecewise linear function (order==2).
1224
1225
1226    Precondition:
1227      All quantities defined at centroids and bed elevation defined at
1228      vertices.
1229
1230    Postcondition
1231      Conserved quantities defined at vertices
1232    """
1233
1234    # Remove very thin layers of water
1235    protect_against_infinitesimal_and_negative_heights(domain)
1236
1237    for name in domain.conserved_quantities:
1238        Q = domain.quantities[name]
1239        if domain._order_ == 1:
1240            Q.extrapolate_first_order()
1241        elif domain._order_ == 2:
1242            Q.extrapolate_second_order_and_limit_by_edge()
1243        else:
1244            raise 'Unknown order'
1245
1246    balance_deep_and_shallow(domain)
1247
1248    # Compute edge values by interpolation
1249    for name in domain.conserved_quantities:
1250        Q = domain.quantities[name]
1251        Q.interpolate_from_vertices_to_edges()
1252
1253##
1254# @brief  Protect against infinitesimal heights and associated high velocities.
1255# @param domain The domain to operate on.
1256def protect_against_infinitesimal_and_negative_heights(domain):
1257    """Protect against infinitesimal heights and associated high velocities"""
1258
1259    from shallow_water_ext import protect
1260
1261    # Shortcuts
1262    wc = domain.quantities['stage'].centroid_values
1263    zc = domain.quantities['elevation'].centroid_values
1264    xmomc = domain.quantities['xmomentum'].centroid_values
1265    ymomc = domain.quantities['ymomentum'].centroid_values
1266
1267    protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
1268            domain.epsilon, wc, zc, xmomc, ymomc)
1269
1270##
1271# @brief Wrapper for C function balance_deep_and_shallow_c().
1272# @param domain The domain to operate on.
1273def balance_deep_and_shallow(domain):
1274    """Compute linear combination between stage as computed by
1275    gradient-limiters limiting using w, and stage computed by
1276    gradient-limiters limiting using h (h-limiter).
1277    The former takes precedence when heights are large compared to the
1278    bed slope while the latter takes precedence when heights are
1279    relatively small.  Anything in between is computed as a balanced
1280    linear combination in order to avoid numerical disturbances which
1281    would otherwise appear as a result of hard switching between
1282    modes.
1283
1284    Wrapper for C implementation
1285    """
1286
1287    from shallow_water_ext import balance_deep_and_shallow \
1288                                  as balance_deep_and_shallow_c
1289
1290    # Shortcuts
1291    wc = domain.quantities['stage'].centroid_values
1292    zc = domain.quantities['elevation'].centroid_values
1293    wv = domain.quantities['stage'].vertex_values
1294    zv = domain.quantities['elevation'].vertex_values
1295
1296    # Momentums at centroids
1297    xmomc = domain.quantities['xmomentum'].centroid_values
1298    ymomc = domain.quantities['ymomentum'].centroid_values
1299
1300    # Momentums at vertices
1301    xmomv = domain.quantities['xmomentum'].vertex_values
1302    ymomv = domain.quantities['ymomentum'].vertex_values
1303
1304    balance_deep_and_shallow_c(domain,
1305                               wc, zc, wv, zv, wc,
1306                               xmomc, ymomc, xmomv, ymomv)
1307
1308
1309################################################################################
1310# Boundary conditions - specific to the shallow water wave equation
1311################################################################################
1312
1313##
1314# @brief Class for a reflective boundary.
1315# @note Inherits from Boundary.
1316class Reflective_boundary(Boundary):
1317    """Reflective boundary returns same conserved quantities as
1318    those present in its neighbour volume but reflected.
1319
1320    This class is specific to the shallow water equation as it
1321    works with the momentum quantities assumed to be the second
1322    and third conserved quantities.
1323    """
1324
1325    ##
1326    # @brief Instantiate a Reflective_boundary.
1327    # @param domain
1328    def __init__(self, domain=None):
1329        Boundary.__init__(self)
1330
1331        if domain is None:
1332            msg = 'Domain must be specified for reflective boundary'
1333            raise Exception, msg
1334
1335        # Handy shorthands
1336        self.stage = domain.quantities['stage'].edge_values
1337        self.xmom = domain.quantities['xmomentum'].edge_values
1338        self.ymom = domain.quantities['ymomentum'].edge_values
1339        self.normals = domain.normals
1340
1341        self.conserved_quantities = num.zeros(3, num.float)
1342
1343    ##
1344    # @brief Return a representation of this instance.
1345    def __repr__(self):
1346        return 'Reflective_boundary'
1347
1348    ##
1349    # @brief Calculate reflections (reverse outward momentum).
1350    # @param vol_id
1351    # @param edge_id
1352    def evaluate(self, vol_id, edge_id):
1353        """Reflective boundaries reverses the outward momentum
1354        of the volume they serve.
1355        """
1356
1357        q = self.conserved_quantities
1358        q[0] = self.stage[vol_id, edge_id]
1359        q[1] = self.xmom[vol_id, edge_id]
1360        q[2] = self.ymom[vol_id, edge_id]
1361
1362        normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
1363
1364        r = rotate(q, normal, direction = 1)
1365        r[1] = -r[1]
1366        q = rotate(r, normal, direction = -1)
1367
1368        return q
1369
1370
1371##
1372# @brief Class for a transmissive boundary.
1373# @note Inherits from Boundary.
1374class Transmissive_momentum_set_stage_boundary(Boundary):
1375    """Returns same momentum conserved quantities as
1376    those present in its neighbour volume.
1377    Sets stage by specifying a function f of time which may either be a
1378    vector function or a scalar function
1379
1380    Example:
1381
1382    def waveform(t):
1383        return sea_level + normalized_amplitude/cosh(t-25)**2
1384
1385    Bts = Transmissive_momentum_set_stage_boundary(domain, waveform)
1386
1387    Underlying domain must be specified when boundary is instantiated
1388    """
1389
1390    ##
1391    # @brief Instantiate a Transmissive_momentum_set_stage_boundary.
1392    # @param domain
1393    # @param function
1394    def __init__(self, domain=None, function=None):
1395        Boundary.__init__(self)
1396
1397        if domain is None:
1398            msg = 'Domain must be specified for this type boundary'
1399            raise Exception, msg
1400
1401        if function is None:
1402            msg = 'Function must be specified for this type boundary'
1403            raise Exception, msg
1404
1405        self.domain = domain
1406        self.function = function
1407
1408    ##
1409    # @brief Return a representation of this instance.
1410    def __repr__(self):
1411        return 'Transmissive_momentum_set_stage_boundary(%s)' %self.domain
1412
1413    ##
1414    # @brief Calculate transmissive results.
1415    # @param vol_id
1416    # @param edge_id
1417    def evaluate(self, vol_id, edge_id):
1418        """Transmissive momentum set stage boundaries return the edge momentum
1419        values of the volume they serve.
1420        """
1421
1422        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
1423
1424        t = self.domain.get_time()
1425
1426        if hasattr(self.function, 'time'):
1427            # Roll boundary over if time exceeds
1428            while t > self.function.time[-1]:
1429                msg = 'WARNING: domain time %.2f has exceeded' % t
1430                msg += 'time provided in '
1431                msg += 'transmissive_momentum_set_stage_boundary object.\n'
1432                msg += 'I will continue, reusing the object from t==0'
1433                log.critical(msg)
1434                t -= self.function.time[-1]
1435
1436        value = self.function(t)
1437        try:
1438            x = float(value)
1439        except:
1440            x = float(value[0])
1441
1442        q[0] = x
1443           
1444        return q
1445
1446        # FIXME: Consider this (taken from File_boundary) to allow
1447        # spatial variation
1448        # if vol_id is not None and edge_id is not None:
1449        #     i = self.boundary_indices[ vol_id, edge_id ]
1450        #     return self.F(t, point_id = i)
1451        # else:
1452        #     return self.F(t)
1453
1454
1455##
1456# @brief Deprecated boundary class.
1457class Transmissive_Momentum_Set_Stage_boundary(Transmissive_momentum_set_stage_boundary):
1458    pass
1459
1460
1461##
1462# @brief Class for a transmissive boundary.
1463# @note Inherits from Boundary.
1464class Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(Boundary):
1465    """Returns the same normal momentum as that
1466    present in neighbour volume edge. Zero out the tangential momentum.
1467    Sets stage by specifying a function f of time which may either be a
1468    vector function or a scalar function
1469
1470    Example:
1471
1472    def waveform(t):
1473        return sea_level + normalized_amplitude/cosh(t-25)**2
1474
1475    Bts = Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(domain, waveform)
1476
1477    Underlying domain must be specified when boundary is instantiated
1478    """
1479
1480    ##
1481    # @brief Instantiate a Transmissive_n_momentum_zero_t_momentum_set_stage_boundary.
1482    # @param domain
1483    # @param function
1484    def __init__(self, domain=None, function=None):
1485        Boundary.__init__(self)
1486
1487        if domain is None:
1488            msg = 'Domain must be specified for this type boundary'
1489            raise Exception, msg
1490
1491        if function is None:
1492            msg = 'Function must be specified for this type boundary'
1493            raise Exception, msg
1494
1495        self.domain = domain
1496        self.function = function
1497
1498    ##
1499    # @brief Return a representation of this instance.
1500    def __repr__(self):
1501        return 'Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(%s)' %self.domain
1502
1503    ##
1504    # @brief Calculate transmissive results.
1505    # @param vol_id
1506    # @param edge_id
1507    def evaluate(self, vol_id, edge_id):
1508        """Transmissive_n_momentum_zero_t_momentum_set_stage_boundary
1509        return the edge momentum values of the volume they serve.
1510        """
1511
1512        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
1513
1514        normal = self.domain.get_normal(vol_id, edge_id)
1515
1516
1517        t = self.domain.get_time()
1518
1519        if hasattr(self.function, 'time'):
1520            # Roll boundary over if time exceeds
1521            while t > self.function.time[-1]:
1522                msg = 'WARNING: domain time %.2f has exceeded' % t
1523                msg += 'time provided in '
1524                msg += 'transmissive_momentum_set_stage_boundary object.\n'
1525                msg += 'I will continue, reusing the object from t==0'
1526                log.critical(msg)
1527                t -= self.function.time[-1]
1528
1529        value = self.function(t)
1530        try:
1531            x = float(value)
1532        except:
1533            x = float(value[0])
1534
1535        ## import math
1536        ## if vol_id == 9433:
1537        ##     print 'vol_id = ',vol_id, ' edge_id = ',edge_id, 'q = ', q, ' x = ',x
1538        ##     print 'normal = ', normal
1539        ##     print 'n . p = ', (normal[0]*q[1] + normal[1]*q[2])
1540        ##     print 't . p = ', (normal[1]*q[1] - normal[0]*q[2])
1541
1542
1543        q[0] = x
1544        ndotq = (normal[0]*q[1] + normal[1]*q[2])
1545        q[1] = normal[0]*ndotq
1546        q[2] = normal[1]*ndotq
1547
1548           
1549        return q
1550
1551##
1552# @brief A transmissive boundary, momentum set to zero.
1553# @note Inherits from Bouondary.
1554class Transmissive_stage_zero_momentum_boundary(Boundary):
1555    """Return same stage as those present in its neighbour volume.
1556    Set momentum to zero.
1557
1558    Underlying domain must be specified when boundary is instantiated
1559    """
1560
1561    ##
1562    # @brief Instantiate a Transmissive (zero momentum) boundary.
1563    # @param domain
1564    def __init__(self, domain=None):
1565        Boundary.__init__(self)
1566
1567        if domain is None:
1568            msg = ('Domain must be specified for '
1569                   'Transmissive_stage_zero_momentum boundary')
1570            raise Exception, msg
1571
1572        self.domain = domain
1573
1574    ##
1575    # @brief Return a representation of this instance.
1576    def __repr__(self):
1577        return 'Transmissive_stage_zero_momentum_boundary(%s)' % self.domain
1578
1579    ##
1580    # @brief Calculate transmissive (zero momentum) results.
1581    # @param vol_id
1582    # @param edge_id
1583    def evaluate(self, vol_id, edge_id):
1584        """Transmissive boundaries return the edge values
1585        of the volume they serve.
1586        """
1587
1588        q = self.domain.get_conserved_quantities(vol_id, edge=edge_id)
1589
1590        q[1] = q[2] = 0.0
1591        return q
1592
1593
1594##
1595# @brief Class for a Dirichlet discharge boundary.
1596# @note Inherits from Boundary.
1597class Dirichlet_discharge_boundary(Boundary):
1598    """
1599    Sets stage (stage0)
1600    Sets momentum (wh0) in the inward normal direction.
1601
1602    Underlying domain must be specified when boundary is instantiated
1603    """
1604
1605    ##
1606    # @brief Instantiate a Dirichlet discharge boundary.
1607    # @param domain
1608    # @param stage0
1609    # @param wh0
1610    def __init__(self, domain=None, stage0=None, wh0=None):
1611        Boundary.__init__(self)
1612
1613        if domain is None:
1614            msg = 'Domain must be specified for this type of boundary'
1615            raise Exception, msg
1616
1617        if stage0 is None:
1618            raise Exception, 'Stage must be specified for this type of boundary'
1619
1620        if wh0 is None:
1621            wh0 = 0.0
1622
1623        self.domain = domain
1624        self.stage0 = stage0
1625        self.wh0 = wh0
1626
1627    ##
1628    # @brief Return a representation of this instance.
1629    def __repr__(self):
1630        return 'Dirichlet_Discharge_boundary(%s)' % self.domain
1631
1632    ##
1633    # @brief Calculate Dirichlet discharge boundary results.
1634    # @param vol_id
1635    # @param edge_id
1636    def evaluate(self, vol_id, edge_id):
1637        """Set discharge in the (inward) normal direction"""
1638
1639        normal = self.domain.get_normal(vol_id,edge_id)
1640        q = [self.stage0, -self.wh0*normal[0], -self.wh0*normal[1]]
1641        return q
1642
1643        # FIXME: Consider this (taken from File_boundary) to allow
1644        # spatial variation
1645        # if vol_id is not None and edge_id is not None:
1646        #     i = self.boundary_indices[ vol_id, edge_id ]
1647        #     return self.F(t, point_id = i)
1648        # else:
1649        #     return self.F(t)
1650
1651
1652# Backward compatibility
1653# FIXME(Ole): Deprecate
1654##
1655# @brief Deprecated
1656class Dirichlet_Discharge_boundary(Dirichlet_discharge_boundary):
1657    pass
1658
1659
1660class Inflow_boundary(Boundary):
1661    """Apply given flow in m^3/s to boundary segment.
1662    Depth and momentum is derived using Manning's formula.
1663
1664    Underlying domain must be specified when boundary is instantiated
1665    """
1666   
1667    # FIXME (Ole): This is work in progress and definitely not finished.
1668    # The associated test has been disabled
1669
1670    def __init__(self, domain=None, rate=0.0):
1671        Boundary.__init__(self)
1672
1673        if domain is None:
1674            msg = 'Domain must be specified for '
1675            msg += 'Inflow boundary'
1676            raise Exception, msg
1677
1678        self.domain = domain
1679       
1680        # FIXME(Ole): Allow rate to be time dependent as well
1681        self.rate = rate
1682        self.tag = None # Placeholder for tag associated with this object.
1683
1684    def __repr__(self):
1685        return 'Inflow_boundary(%s)' %self.domain
1686
1687    def evaluate(self, vol_id, edge_id):
1688        """Apply inflow rate at each edge of this boundary
1689        """
1690       
1691        # First find all segments having the same tag is vol_id, edge_id
1692        # This will be done the first time evaluate is called.
1693        if self.tag is None:
1694            boundary = self.domain.boundary
1695            self.tag = boundary[(vol_id, edge_id)]       
1696           
1697            # Find total length of boundary with this tag
1698            length = 0.0
1699            for v_id, e_id in boundary:
1700                if self.tag == boundary[(v_id, e_id)]:
1701                    length += self.domain.mesh.get_edgelength(v_id, e_id)           
1702
1703            self.length = length
1704            self.average_momentum = self.rate/length
1705           
1706           
1707        # Average momentum has now been established across this boundary
1708        # Compute momentum in the inward normal direction
1709       
1710        inward_normal = -self.domain.mesh.get_normal(vol_id, edge_id)       
1711        xmomentum, ymomentum = self.average_momentum * inward_normal
1712           
1713        # Compute depth based on Manning's formula v = 1/n h^{2/3} sqrt(S)
1714        # Where v is velocity, n is manning's coefficient, h is depth and S is the slope into the domain.
1715        # Let mu be the momentum (vh), then this equation becomes: mu = 1/n h^{5/3} sqrt(S)
1716        # from which we can isolate depth to get
1717        # h = (mu n/sqrt(S) )^{3/5}
1718       
1719        slope = 0 # get gradient for this triangle dot normal
1720       
1721        # get manning coef from this triangle
1722        friction = self.domain.get_quantity('friction').get_values(location='edges', 
1723                                                                   indices=[vol_id])[0]
1724        mannings_n = friction[edge_id]
1725
1726        if slope > epsilon and mannings_n > epsilon:
1727            depth = pow(self.average_momentum * mannings_n/math.sqrt(slope), 3.0/5) 
1728        else:
1729            depth = 1.0
1730           
1731        # Elevation on this edge   
1732       
1733        z = self.domain.get_quantity('elevation').get_values(location='edges', 
1734                                                             indices=[vol_id])[0]
1735        elevation = z[edge_id]
1736           
1737        # Assign conserved quantities and return
1738        q = num.array([elevation + depth, xmomentum, ymomentum], num.Float)
1739        return q
1740
1741
1742       
1743   
1744           
1745       
1746class Field_boundary(Boundary):
1747    """Set boundary from given field represented in an sww file containing
1748    values for stage, xmomentum and ymomentum.
1749
1750    Optionally, the user can specify mean_stage to offset the stage provided
1751    in the sww file.
1752
1753    This function is a thin wrapper around the generic File_boundary. The
1754    difference between the file_boundary and field_boundary is only that the
1755    field_boundary will allow you to change the level of the stage height when
1756    you read in the boundary condition. This is very useful when running
1757    different tide heights in the same area as you need only to convert one
1758    boundary condition to a SWW file, ideally for tide height of 0 m
1759    (saving disk space). Then you can use field_boundary to read this SWW file
1760    and change the stage height (tide) on the fly depending on the scenario.
1761    """
1762
1763    ##
1764    # @brief Construct an instance of a 'field' boundary.
1765    # @param filename Name of SWW file containing stage, x and ymomentum.
1766    # @param domain Shallow water domain for which the boundary applies.
1767    # @param mean_stage Mean water level added to stage derived from SWW file.
1768    # @param time_thinning Time step thinning factor.
1769    # @param time_limit
1770    # @param boundary_polygon
1771    # @param default_boundary None or an instance of Boundary.
1772    # @param use_cache True if caching is to be used.
1773    # @param verbose True if this method is to be verbose.
1774    def __init__(self,
1775                 filename,
1776                 domain,
1777                 mean_stage=0.0,
1778                 time_thinning=1,
1779                 time_limit=None,
1780                 boundary_polygon=None,
1781                 default_boundary=None,
1782                 use_cache=False,
1783                 verbose=False):
1784        """Constructor
1785
1786        filename: Name of sww file
1787        domain: pointer to shallow water domain for which the boundary applies
1788        mean_stage: The mean water level which will be added to stage derived
1789                    from the boundary condition
1790        time_thinning: Will set how many time steps from the sww file read in
1791                       will be interpolated to the boundary. For example if
1792                       the sww file has 1 second time steps and is 24 hours
1793                       in length it has 86400 time steps. If you set
1794                       time_thinning to 1 it will read all these steps.
1795                       If you set it to 100 it will read every 100th step eg
1796                       only 864 step. This parameter is very useful to increase
1797                       the speed of a model run that you are setting up
1798                       and testing.
1799
1800        default_boundary: Must be either None or an instance of a
1801                          class descending from class Boundary.
1802                          This will be used in case model time exceeds
1803                          that available in the underlying data.
1804
1805                          Note that mean_stage will also be added to this.
1806                                               
1807        use_cache:
1808        verbose:
1809        """
1810
1811        # Create generic file_boundary object
1812        self.file_boundary = File_boundary(filename,
1813                                           domain,
1814                                           time_thinning=time_thinning,
1815                                           time_limit=time_limit,
1816                                           boundary_polygon=boundary_polygon,
1817                                           default_boundary=default_boundary,
1818                                           use_cache=use_cache,
1819                                           verbose=verbose)
1820
1821        # Record information from File_boundary
1822        self.F = self.file_boundary.F
1823        self.domain = self.file_boundary.domain
1824
1825        # Record mean stage
1826        self.mean_stage = mean_stage
1827
1828    ##
1829    # @note Generate a string representation of this instance.
1830    def __repr__(self):
1831        return 'Field boundary'
1832
1833    ##
1834    # @brief Calculate 'field' boundary results.
1835    # @param vol_id
1836    # @param edge_id
1837    def evaluate(self, vol_id=None, edge_id=None):
1838        """Return linearly interpolated values based on domain.time
1839
1840        vol_id and edge_id are ignored
1841        """
1842
1843        # Evaluate file boundary
1844        q = self.file_boundary.evaluate(vol_id, edge_id)
1845
1846        # Adjust stage
1847        for j, name in enumerate(self.domain.conserved_quantities):
1848            if name == 'stage':
1849                q[j] += self.mean_stage
1850        return q
1851
1852
1853################################################################################
1854# Standard forcing terms
1855################################################################################
1856
1857##
1858# @brief Apply gravitational pull in the presence of bed slope.
1859# @param domain The domain to apply gravity to.
1860# @note Wrapper for C function gravity_c().
1861def gravity(domain):
1862    """Apply gravitational pull in the presence of bed slope
1863    Wrapper calls underlying C implementation
1864    """
1865
1866    from shallow_water_ext import gravity as gravity_c
1867
1868    xmom = domain.quantities['xmomentum'].explicit_update
1869    ymom = domain.quantities['ymomentum'].explicit_update
1870
1871    stage = domain.quantities['stage']
1872    elevation = domain.quantities['elevation']
1873
1874    h = stage.centroid_values - elevation.centroid_values
1875    z = elevation.vertex_values
1876
1877    x = domain.get_vertex_coordinates()
1878    g = domain.g
1879
1880    gravity_c(g, h, z, x, xmom, ymom)    #, 1.0e-6)
1881
1882##
1883# @brief Apply friction to a surface (implicit).
1884# @param domain The domain to apply Manning friction to.
1885# @note Wrapper for C function manning_friction_c().
1886def manning_friction_implicit(domain):
1887    """Apply (Manning) friction to water momentum
1888    Wrapper for c version
1889    """
1890
1891    from shallow_water_ext import manning_friction as manning_friction_c
1892
1893    xmom = domain.quantities['xmomentum']
1894    ymom = domain.quantities['ymomentum']
1895
1896    w = domain.quantities['stage'].centroid_values
1897    z = domain.quantities['elevation'].centroid_values
1898
1899    uh = xmom.centroid_values
1900    vh = ymom.centroid_values
1901    eta = domain.quantities['friction'].centroid_values
1902
1903    xmom_update = xmom.semi_implicit_update
1904    ymom_update = ymom.semi_implicit_update
1905
1906    N = len(domain)
1907    eps = domain.minimum_allowed_height
1908    g = domain.g
1909
1910    manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
1911
1912##
1913# @brief Apply friction to a surface (explicit).
1914# @param domain The domain to apply Manning friction to.
1915# @note Wrapper for C function manning_friction_c().
1916def manning_friction_explicit(domain):
1917    """Apply (Manning) friction to water momentum
1918    Wrapper for c version
1919    """
1920
1921    from shallow_water_ext import manning_friction as manning_friction_c
1922
1923    xmom = domain.quantities['xmomentum']
1924    ymom = domain.quantities['ymomentum']
1925
1926    w = domain.quantities['stage'].centroid_values
1927    z = domain.quantities['elevation'].centroid_values
1928
1929    uh = xmom.centroid_values
1930    vh = ymom.centroid_values
1931    eta = domain.quantities['friction'].centroid_values
1932
1933    xmom_update = xmom.explicit_update
1934    ymom_update = ymom.explicit_update
1935
1936    N = len(domain)
1937    eps = domain.minimum_allowed_height
1938    g = domain.g
1939
1940    manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
1941
1942
1943# FIXME (Ole): This was implemented for use with one of the analytical solutions (Sampson?)
1944##
1945# @brief Apply linear friction to a surface.
1946# @param domain The domain to apply Manning friction to.
1947# @note Is this still used (30 Oct 2007)?
1948def linear_friction(domain):
1949    """Apply linear friction to water momentum
1950
1951    Assumes quantity: 'linear_friction' to be present
1952    """
1953
1954    from math import sqrt
1955
1956    w = domain.quantities['stage'].centroid_values
1957    z = domain.quantities['elevation'].centroid_values
1958    h = w-z
1959
1960    uh = domain.quantities['xmomentum'].centroid_values
1961    vh = domain.quantities['ymomentum'].centroid_values
1962    tau = domain.quantities['linear_friction'].centroid_values
1963
1964    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1965    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1966
1967    N = len(domain) # number_of_triangles
1968    eps = domain.minimum_allowed_height
1969    g = domain.g #Not necessary? Why was this added?
1970
1971    for k in range(N):
1972        if tau[k] >= eps:
1973            if h[k] >= eps:
1974                S = -tau[k]/h[k]
1975
1976                #Update momentum
1977                xmom_update[k] += S*uh[k]
1978                ymom_update[k] += S*vh[k]
1979
1980def depth_dependent_friction(domain, default_friction,
1981                             surface_roughness_data,
1982                             verbose=False):
1983    """Returns an array of friction values for each wet element adjusted for depth.
1984
1985    Inputs:
1986        domain - computational domain object
1987        default_friction - depth independent bottom friction
1988        surface_roughness_data - N x 5 array of n0, d1, n1, d2, n2 values for each
1989        friction region.
1990
1991    Outputs:
1992        wet_friction - Array that can be used directly to update friction as follows:
1993                       domain.set_quantity('friction', wet_friction)
1994
1995       
1996       
1997    """
1998
1999    import numpy as num
2000   
2001    # Create a temp array to store updated depth dependent friction for wet elements
2002    # EHR this is outwardly inneficient but not obvious how to avoid recreating each call??????
2003    N=len(domain)
2004    wet_friction    = num.zeros(N, num.float)
2005    wet_friction[:] = default_n0   # Initially assign default_n0 to all array so sure have no zeros values
2006   
2007   
2008    depth = domain.create_quantity_from_expression('stage - elevation')  # create depth instance for this timestep
2009    # Recompute depth as vector 
2010    d = depth.get_values(location='centroids')
2011 
2012    # rebuild the 'friction' values adjusted for depth at this instant
2013    for i in domain.get_wet_elements():                                  # loop for each wet element in domain
2014       
2015        # Get roughness data for each element
2016        n0 = float(surface_roughness_data[i,0])
2017        d1 = float(surface_roughness_data[i,1])
2018        n1 = float(surface_roughness_data[i,2])
2019        d2 = float(surface_roughness_data[i,3])
2020        n2 = float(surface_roughness_data[i,4])
2021       
2022       
2023        # Recompute friction values from depth for this element
2024               
2025        if d[i]   <= d1: 
2026            depth_dependent_friction = n1
2027        elif d[i] >= d2:
2028            depth_dependent_friction = n2
2029        else:
2030            depth_dependent_friction = n1+((n2-n1)/(d2-d1))*(d[i]-d1)
2031           
2032        # check sanity of result
2033        if (depth_dependent_friction  < 0.010 or depth_dependent_friction > 9999.0) :
2034            log.critical('%s >>>> WARNING: computed depth_dependent friction '
2035                         'out of range, ddf%f, n1=%f, n2=%f'
2036                         % (model_data.basename,
2037                            depth_dependent_friction, n1, n2))
2038       
2039        # update depth dependent friction  for that wet element
2040        wet_friction[i] = depth_dependent_friction
2041       
2042    # EHR add code to show range of 'friction across domain at this instant as sanity check?????????
2043   
2044    if verbose :
2045        nvals=domain.get_quantity('friction').get_values(location='centroids')        # return array of domain nvals
2046        n_min=min(nvals)
2047        n_max=max(nvals)
2048       
2049        log.critical('         ++++ calculate_depth_dependent_friction - '
2050                     'Updated friction - range  %7.3f to %7.3f'
2051                     % (n_min, n_max))
2052   
2053    return wet_friction
2054
2055
2056################################################################################
2057# Experimental auxiliary functions
2058################################################################################
2059
2060##
2061# @brief Check forcefield parameter.
2062# @param f Object to check.
2063# @note 'f' may be a callable object or a scalar value.
2064def check_forcefield(f):
2065    """Check that force object is as expected.
2066   
2067    Check that f is either:
2068    1: a callable object f(t,x,y), where x and y are vectors
2069       and that it returns an array or a list of same length
2070       as x and y
2071    2: a scalar
2072    """
2073
2074    if callable(f):
2075        N = 3
2076        x = num.ones(3, num.float)
2077        y = num.ones(3, num.float)
2078        try:
2079            q = f(1.0, x=x, y=y)
2080        except Exception, e:
2081            msg = 'Function %s could not be executed:\n%s' %(f, e)
2082            # FIXME: Reconsider this semantics
2083            raise Exception, msg
2084
2085        try:
2086            q = num.array(q, num.float)
2087        except:
2088            msg = ('Return value from vector function %s could not '
2089                   'be converted into a numeric array of floats.\nSpecified '
2090                   'function should return either list or array.' % f)
2091            raise Exception, msg
2092
2093        # Is this really what we want?
2094        # info is "(func name, filename, defining line)"
2095        func_info = (f.func_name, f.func_code.co_filename,
2096                     f.func_code.co_firstlineno)
2097        func_msg = 'Function %s (defined in %s, line %d)' % func_info
2098        try:
2099            result_len = len(q)
2100        except:
2101            msg = '%s must return vector' % func_msg
2102            self.fail(msg)
2103        msg = '%s must return vector of length %d' % (func_msg, N)
2104        assert result_len == N, msg
2105    else:
2106        try:
2107            f = float(f)
2108        except:
2109            msg = ('Force field %s must be a scalar value coercible to float.'
2110                   % str(f))
2111            raise Exception, msg
2112
2113    return f
2114
2115
2116##
2117# Class to apply a wind stress to a domain.
2118class Wind_stress:
2119    """Apply wind stress to water momentum in terms of
2120    wind speed [m/s] and wind direction [degrees]
2121    """
2122
2123    ##
2124    # @brief Create an instance of Wind_stress.
2125    # @param *args
2126    # @param **kwargs
2127    def __init__(self, *args, **kwargs):
2128        """Initialise windfield from wind speed s [m/s]
2129        and wind direction phi [degrees]
2130
2131        Inputs v and phi can be either scalars or Python functions, e.g.
2132
2133        W = Wind_stress(10, 178)
2134
2135        #FIXME - 'normal' degrees are assumed for now, i.e. the
2136        vector (1,0) has zero degrees.
2137        We may need to convert from 'compass' degrees later on and also
2138        map from True north to grid north.
2139
2140        Arguments can also be Python functions of t,x,y as in
2141
2142        def speed(t,x,y):
2143            ...
2144            return s
2145
2146        def angle(t,x,y):
2147            ...
2148            return phi
2149
2150        where x and y are vectors.
2151
2152        and then pass the functions in
2153
2154        W = Wind_stress(speed, angle)
2155
2156        The instantiated object W can be appended to the list of
2157        forcing_terms as in
2158
2159        Alternatively, one vector valued function for (speed, angle)
2160        can be applied, providing both quantities simultaneously.
2161        As in
2162        W = Wind_stress(F), where returns (speed, angle) for each t.
2163
2164        domain.forcing_terms.append(W)
2165        """
2166
2167        from anuga.config import rho_a, rho_w, eta_w
2168
2169        if len(args) == 2:
2170            s = args[0]
2171            phi = args[1]
2172        elif len(args) == 1:
2173            # Assume vector function returning (s, phi)(t,x,y)
2174            vector_function = args[0]
2175            s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
2176            phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
2177        else:
2178           # Assume info is in 2 keyword arguments
2179           if len(kwargs) == 2:
2180               s = kwargs['s']
2181               phi = kwargs['phi']
2182           else:
2183               raise Exception, 'Assumes two keyword arguments: s=..., phi=....'
2184
2185        self.speed = check_forcefield(s)
2186        self.phi = check_forcefield(phi)
2187
2188        self.const = eta_w*rho_a/rho_w
2189
2190    ##
2191    # @brief 'execute' this class instance.
2192    # @param domain
2193    def __call__(self, domain):
2194        """Evaluate windfield based on values found in domain"""
2195
2196        from math import pi, cos, sin, sqrt
2197
2198        xmom_update = domain.quantities['xmomentum'].explicit_update
2199        ymom_update = domain.quantities['ymomentum'].explicit_update
2200
2201        N = len(domain)    # number_of_triangles
2202        t = domain.time
2203
2204        if callable(self.speed):
2205            xc = domain.get_centroid_coordinates()
2206            s_vec = self.speed(t, xc[:,0], xc[:,1])
2207        else:
2208            # Assume s is a scalar
2209            try:
2210                s_vec = self.speed * num.ones(N, num.float)
2211            except:
2212                msg = 'Speed must be either callable or a scalar: %s' %self.s
2213                raise msg
2214
2215        if callable(self.phi):
2216            xc = domain.get_centroid_coordinates()
2217            phi_vec = self.phi(t, xc[:,0], xc[:,1])
2218        else:
2219            # Assume phi is a scalar
2220
2221            try:
2222                phi_vec = self.phi * num.ones(N, num.float)
2223            except:
2224                msg = 'Angle must be either callable or a scalar: %s' %self.phi
2225                raise msg
2226
2227        assign_windfield_values(xmom_update, ymom_update,
2228                                s_vec, phi_vec, self.const)
2229
2230
2231##
2232# @brief Assign wind field values
2233# @param xmom_update
2234# @param ymom_update
2235# @param s_vec
2236# @param phi_vec
2237# @param const
2238def assign_windfield_values(xmom_update, ymom_update,
2239                            s_vec, phi_vec, const):
2240    """Python version of assigning wind field to update vectors.
2241    A C version also exists (for speed)
2242    """
2243
2244    from math import pi, cos, sin, sqrt
2245
2246    N = len(s_vec)
2247    for k in range(N):
2248        s = s_vec[k]
2249        phi = phi_vec[k]
2250
2251        # Convert to radians
2252        phi = phi*pi/180
2253
2254        # Compute velocity vector (u, v)
2255        u = s*cos(phi)
2256        v = s*sin(phi)
2257
2258        # Compute wind stress
2259        S = const * sqrt(u**2 + v**2)
2260        xmom_update[k] += S*u
2261        ymom_update[k] += S*v
2262
2263
2264##
2265# @brief A class for a general explicit forcing term.
2266class General_forcing:
2267    """General explicit forcing term for update of quantity
2268
2269    This is used by Inflow and Rainfall for instance
2270
2271
2272    General_forcing(quantity_name, rate, center, radius, polygon)
2273
2274    domain:     ANUGA computational domain
2275    quantity_name: Name of quantity to update.
2276                   It must be a known conserved quantity.
2277
2278    rate [?/s]: Total rate of change over the specified area.
2279                This parameter can be either a constant or a
2280                function of time. Positive values indicate increases,
2281                negative values indicate decreases.
2282                Rate can be None at initialisation but must be specified
2283                before forcing term is applied (i.e. simulation has started).
2284
2285    center [m]: Coordinates at center of flow point
2286    radius [m]: Size of circular area
2287    polygon:    Arbitrary polygon
2288    default_rate: Rate to be used if rate fails (e.g. if model time exceeds its data)
2289                  Admissible types: None, constant number or function of t
2290
2291
2292    Either center, radius or polygon can be specified but not both.
2293    If neither are specified the entire domain gets updated.
2294    All coordinates to be specified in absolute UTM coordinates (x, y) assuming the zone of domain.
2295
2296    Inflow or Rainfall for examples of use
2297    """
2298
2299
2300    # FIXME (AnyOne) : Add various methods to allow spatial variations
2301
2302    ##
2303    # @brief Create an instance of this forcing term.
2304    # @param domain
2305    # @param quantity_name
2306    # @param rate
2307    # @param center
2308    # @param radius
2309    # @param polygon
2310    # @param default_rate
2311    # @param verbose
2312    def __init__(self,
2313                 domain,
2314                 quantity_name,
2315                 rate=0.0,
2316                 center=None,
2317                 radius=None,
2318                 polygon=None,
2319                 default_rate=None,
2320                 verbose=False):
2321
2322        from math import pi, cos, sin
2323
2324        if center is None:
2325            msg = 'I got radius but no center.'
2326            assert radius is None, msg
2327
2328        if radius is None:
2329            msg += 'I got center but no radius.'
2330            assert center is None, msg
2331
2332        self.domain = domain
2333        self.quantity_name = quantity_name
2334        self.rate = rate
2335        self.center = ensure_numeric(center)
2336        self.radius = radius
2337        self.polygon = polygon
2338        self.verbose = verbose
2339        self.value = 0.0    # Can be used to remember value at
2340                            # previous timestep in order to obtain rate
2341
2342        # Get boundary (in absolute coordinates)
2343        bounding_polygon = domain.get_boundary_polygon()
2344
2345        # Update area if applicable
2346        if center is not None and radius is not None:
2347            assert len(center) == 2
2348            msg = 'Polygon cannot be specified when center and radius are'
2349            assert polygon is None, msg
2350
2351            # Check that circle center lies within the mesh.
2352            msg = 'Center %s specified for forcing term did not' % str(center)
2353            msg += 'fall within the domain boundary.'
2354            assert is_inside_polygon(center, bounding_polygon), msg
2355
2356            # Check that circle periphery lies within the mesh.
2357            N = 100
2358            periphery_points = []
2359            for i in range(N):
2360                theta = 2*pi*i/100
2361
2362                x = center[0] + radius*cos(theta)
2363                y = center[1] + radius*sin(theta)
2364
2365                periphery_points.append([x,y])
2366
2367            for point in periphery_points:
2368                msg = 'Point %s on periphery for forcing term' % str(point)
2369                msg += ' did not fall within the domain boundary.'
2370                assert is_inside_polygon(point, bounding_polygon), msg
2371
2372        if polygon is not None:
2373            # Check that polygon lies within the mesh.
2374            for point in self.polygon:
2375                msg = 'Point %s in polygon for forcing term' % str(point)
2376                msg += ' did not fall within the domain boundary.'
2377                assert is_inside_polygon(point, bounding_polygon), msg
2378
2379        # Pointer to update vector
2380        self.update = domain.quantities[self.quantity_name].explicit_update
2381
2382        # Determine indices in flow area
2383        N = len(domain)
2384        points = domain.get_centroid_coordinates(absolute=True)
2385
2386        # Calculate indices in exchange area for this forcing term
2387        self.exchange_indices = None
2388        if self.center is not None and self.radius is not None:
2389            # Inlet is circular
2390            inlet_region = 'center=%s, radius=%s' % (self.center, self.radius)
2391
2392            self.exchange_indices = []
2393            for k in range(N):
2394                x, y = points[k,:]    # Centroid
2395
2396                c = self.center
2397                if ((x-c[0])**2+(y-c[1])**2) < self.radius**2:
2398                    self.exchange_indices.append(k)
2399
2400        if self.polygon is not None:
2401            # Inlet is polygon
2402            inlet_region = 'polygon=%s' % (self.polygon) 
2403            self.exchange_indices = inside_polygon(points, self.polygon)
2404
2405        if self.exchange_indices is None:
2406            self.exchange_area = polygon_area(bounding_polygon)
2407        else:   
2408            if len(self.exchange_indices) == 0:
2409                msg = 'No triangles have been identified in '
2410                msg += 'specified region: %s' % inlet_region
2411                raise Exception, msg
2412
2413            # Compute exchange area as the sum of areas of triangles identified
2414            # by circle or polygon
2415            self.exchange_area = 0.0
2416            for i in self.exchange_indices:
2417                self.exchange_area += domain.areas[i]
2418           
2419
2420        msg = 'Exchange area in forcing term'
2421        msg += ' has area = %f' %self.exchange_area
2422        assert self.exchange_area > 0.0           
2423           
2424               
2425
2426           
2427        # Check and store default_rate
2428        msg = ('Keyword argument default_rate must be either None '
2429               'or a function of time.\nI got %s.' % str(default_rate))
2430        assert (default_rate is None or
2431                type(default_rate) in [IntType, FloatType] or
2432                callable(default_rate)), msg
2433
2434        if default_rate is not None:
2435            # If it is a constant, make it a function
2436            if not callable(default_rate):
2437                tmp = default_rate
2438                default_rate = lambda t: tmp
2439
2440            # Check that default_rate is a function of one argument
2441            try:
2442                default_rate(0.0)
2443            except:
2444                raise Exception, msg
2445
2446        self.default_rate = default_rate
2447        self.default_rate_invoked = False    # Flag
2448
2449    ##
2450    # @brief Execute this instance.
2451    # @param domain
2452    def __call__(self, domain):
2453        """Apply inflow function at time specified in domain, update stage"""
2454
2455        # Call virtual method allowing local modifications
2456        t = domain.get_time()
2457        try:
2458            rate = self.update_rate(t)
2459        except Modeltime_too_early, e:
2460            raise Modeltime_too_early, e
2461        except Modeltime_too_late, e:
2462            if self.default_rate is None:
2463                raise Exception, e    # Reraise exception
2464            else:
2465                # Pass control to default rate function
2466                rate = self.default_rate(t)
2467
2468                if self.default_rate_invoked is False:
2469                    # Issue warning the first time
2470                    msg = ('%s\n'
2471                           'Instead I will use the default rate: %s\n'
2472                           'Note: Further warnings will be supressed'
2473                           % (str(e), str(self.default_rate)))
2474                    warn(msg)
2475
2476                    # FIXME (Ole): Replace this crude flag with
2477                    # Python's ability to print warnings only once.
2478                    # See http://docs.python.org/lib/warning-filter.html
2479                    self.default_rate_invoked = True
2480
2481        if rate is None:
2482            msg = ('Attribute rate must be specified in General_forcing '
2483                   'or its descendants before attempting to call it')
2484            raise Exception, msg
2485
2486        # Now rate is a number
2487        if self.verbose is True:
2488            log.critical('Rate of %s at time = %.2f = %f'
2489                         % (self.quantity_name, domain.get_time(), rate))
2490
2491        if self.exchange_indices is None:
2492            self.update[:] += rate
2493        else:
2494            # Brute force assignment of restricted rate
2495            for k in self.exchange_indices:
2496                self.update[k] += rate
2497
2498    ##
2499    # @brief Update the internal rate.
2500    # @param t A callable or scalar used to set the rate.
2501    # @return The new rate.
2502    def update_rate(self, t):
2503        """Virtual method allowing local modifications by writing an
2504        overriding version in descendant
2505        """
2506
2507        if callable(self.rate):
2508            rate = self.rate(t)
2509        else:
2510            rate = self.rate
2511
2512        return rate
2513
2514    ##
2515    # @brief Get values for the specified quantity.
2516    # @param quantity_name Name of the quantity of interest.
2517    # @return The value(s) of the quantity.
2518    # @note If 'quantity_name' is None, use self.quantity_name.
2519    def get_quantity_values(self, quantity_name=None):
2520        """Return values for specified quantity restricted to opening
2521
2522        Optionally a quantity name can be specified if values from another
2523        quantity is sought
2524        """
2525
2526        if quantity_name is None:
2527            quantity_name = self.quantity_name
2528
2529        q = self.domain.quantities[quantity_name]
2530        return q.get_values(location='centroids',
2531                            indices=self.exchange_indices)
2532
2533    ##
2534    # @brief Set value for the specified quantity.
2535    # @param val The value object used to set value.
2536    # @param quantity_name Name of the quantity of interest.
2537    # @note If 'quantity_name' is None, use self.quantity_name.
2538    def set_quantity_values(self, val, quantity_name=None):
2539        """Set values for specified quantity restricted to opening
2540
2541        Optionally a quantity name can be specified if values from another
2542        quantity is sought
2543        """
2544
2545        if quantity_name is None:
2546            quantity_name = self.quantity_name
2547
2548        q = self.domain.quantities[self.quantity_name]
2549        q.set_values(val,
2550                     location='centroids',
2551                     indices=self.exchange_indices)
2552
2553
2554##
2555# @brief A class for rainfall forcing function.
2556# @note Inherits from General_forcing.
2557class Rainfall(General_forcing):
2558    """Class Rainfall - general 'rain over entire domain' forcing term.
2559
2560    Used for implementing Rainfall over the entire domain.
2561
2562        Current Limited to only One Gauge..
2563
2564        Need to add Spatial Varying Capability
2565        (This module came from copying and amending the Inflow Code)
2566
2567    Rainfall(rain)
2568
2569    domain
2570    rain [mm/s]:  Total rain rate over the specified domain.
2571                  NOTE: Raingauge Data needs to reflect the time step.
2572                  IE: if Gauge is mm read at a time step, then the input
2573                  here is as mm/(timeStep) so 10mm in 5minutes becomes
2574                  10/(5x60) = 0.0333mm/s.
2575
2576                  This parameter can be either a constant or a
2577                  function of time. Positive values indicate inflow,
2578                  negative values indicate outflow.
2579                  (and be used for Infiltration - Write Seperate Module)
2580                  The specified flow will be divided by the area of
2581                  the inflow region and then applied to update the
2582                  stage quantity.
2583
2584    polygon: Specifies a polygon to restrict the rainfall.
2585
2586    Examples
2587    How to put them in a run File...
2588
2589    #------------------------------------------------------------------------
2590    # Setup specialised forcing terms
2591    #------------------------------------------------------------------------
2592    # This is the new element implemented by Ole and Rudy to allow direct
2593    # input of Rainfall in mm/s
2594
2595    catchmentrainfall = Rainfall(rain=file_function('Q100_2hr_Rain.tms'))
2596                        # Note need path to File in String.
2597                        # Else assumed in same directory
2598
2599    domain.forcing_terms.append(catchmentrainfall)
2600    """
2601
2602    ##
2603    # @brief Create an instance of the class.
2604    # @param domain Domain of interest.
2605    # @param rate Total rain rate over the specified domain (mm/s).
2606    # @param center
2607    # @param radius
2608    # @param polygon Polygon  to restrict rainfall.
2609    # @param default_rate
2610    # @param verbose True if this instance is to be verbose.
2611    def __init__(self,
2612                 domain,
2613                 rate=0.0,
2614                 center=None,
2615                 radius=None,
2616                 polygon=None,
2617                 default_rate=None,
2618                 verbose=False):
2619
2620        # Converting mm/s to m/s to apply in ANUGA)
2621        if callable(rate):
2622            rain = lambda t: rate(t)/1000.0
2623        else:
2624            rain = rate/1000.0
2625
2626        if default_rate is not None:
2627            if callable(default_rate):
2628                default_rain = lambda t: default_rate(t)/1000.0
2629            else:
2630                default_rain = default_rate/1000.0
2631        else:
2632            default_rain = None
2633
2634           
2635           
2636        General_forcing.__init__(self,
2637                                 domain,
2638                                 'stage',
2639                                 rate=rain,
2640                                 center=center,
2641                                 radius=radius,
2642                                 polygon=polygon,
2643                                 default_rate=default_rain,
2644                                 verbose=verbose)
2645
2646
2647##
2648# @brief A class for inflow (rain and drain) forcing function.
2649# @note Inherits from General_forcing.
2650class Inflow(General_forcing):
2651    """Class Inflow - general 'rain and drain' forcing term.
2652
2653    Useful for implementing flows in and out of the domain.
2654
2655    Inflow(flow, center, radius, polygon)
2656
2657    domain
2658    rate [m^3/s]: Total flow rate over the specified area.
2659                  This parameter can be either a constant or a
2660                  function of time. Positive values indicate inflow,
2661                  negative values indicate outflow.
2662                  The specified flow will be divided by the area of
2663                  the inflow region and then applied to update stage.
2664    center [m]: Coordinates at center of flow point
2665    radius [m]: Size of circular area
2666    polygon:    Arbitrary polygon.
2667
2668    Either center, radius or polygon must be specified
2669
2670    Examples
2671
2672    # Constant drain at 0.003 m^3/s.
2673    # The outflow area is 0.07**2*pi=0.0154 m^2
2674    # This corresponds to a rate of change of 0.003/0.0154 = 0.2 m/s
2675    #
2676    Inflow((0.7, 0.4), 0.07, -0.003)
2677
2678
2679    # Tap turning up to a maximum inflow of 0.0142 m^3/s.
2680    # The inflow area is 0.03**2*pi = 0.00283 m^2
2681    # This corresponds to a rate of change of 0.0142/0.00283 = 5 m/s
2682    # over the specified area
2683    Inflow((0.5, 0.5), 0.03, lambda t: min(0.01*t, 0.0142))
2684
2685
2686    #------------------------------------------------------------------------
2687    # Setup specialised forcing terms
2688    #------------------------------------------------------------------------
2689    # This is the new element implemented by Ole to allow direct input
2690    # of Inflow in m^3/s
2691
2692    hydrograph = Inflow(center=(320, 300), radius=10,
2693                        rate=file_function('Q/QPMF_Rot_Sub13.tms'))
2694
2695    domain.forcing_terms.append(hydrograph)
2696    """
2697
2698    ##
2699    # @brief Create an instance of the class.
2700    # @param domain Domain of interest.
2701    # @param rate Total rain rate over the specified domain (mm/s).
2702    # @param center
2703    # @param radius
2704    # @param polygon Polygon  to restrict rainfall.
2705    # @param default_rate
2706    # @param verbose True if this instance is to be verbose.
2707    def __init__(self,
2708                 domain,
2709                 rate=0.0,
2710                 center=None,
2711                 radius=None,
2712                 polygon=None,
2713                 default_rate=None,
2714                 verbose=False):
2715        # Create object first to make area is available
2716        General_forcing.__init__(self,
2717                                 domain,
2718                                 'stage',
2719                                 rate=rate,
2720                                 center=center,
2721                                 radius=radius,
2722                                 polygon=polygon,
2723                                 default_rate=default_rate,
2724                                 verbose=verbose)
2725
2726    ##
2727    # @brief Update the instance rate.
2728    # @param t New rate object.
2729    def update_rate(self, t):
2730        """Virtual method allowing local modifications by writing an
2731        overriding version in descendant
2732
2733        This one converts m^3/s to m/s which can be added directly
2734        to 'stage' in ANUGA
2735        """
2736
2737        if callable(self.rate):
2738            _rate = self.rate(t)/self.exchange_area
2739        else:
2740            _rate = self.rate/self.exchange_area
2741
2742        return _rate
2743
2744
2745##
2746# @brief A class for creating cross sections.
2747# @note Inherits from General_forcing.
2748class Cross_section:
2749    """Class Cross_section - a class to setup a cross section from
2750    which you can then calculate flow and energy through cross section
2751
2752
2753    Cross_section(domain, polyline)
2754
2755    domain:
2756    polyline: Representation of desired cross section - it may contain
2757              multiple sections allowing for complex shapes. Assume
2758              absolute UTM coordinates.
2759              Format [[x0, y0], [x1, y1], ...]
2760    verbose:
2761    """
2762
2763    ##
2764    # @brief Create an instance of the class.
2765    # @param domain Domain of interest.
2766    # @param polyline Polyline defining cross section
2767    # @param verbose True if this instance is to be verbose.
2768    def __init__(self,
2769                 domain,
2770                 polyline=None,
2771                 verbose=False):
2772       
2773        self.domain = domain
2774        self.polyline = polyline
2775        self.verbose = verbose
2776       
2777        # Find all intersections and associated triangles.
2778        self.segments = self.domain.get_intersecting_segments(self.polyline,
2779                                                              use_cache=True,
2780                                                              verbose=self.verbose)
2781       
2782        # Get midpoints
2783        self.midpoints = segment_midpoints(self.segments)
2784
2785        # Make midpoints Geospatial instances
2786        self.midpoints = ensure_geospatial(self.midpoints, self.domain.geo_reference)
2787
2788
2789    ##
2790    # @brief calculate current flow through cross section
2791    def get_flow_through_cross_section(self):
2792        """ Output: Total flow [m^3/s] across cross section.
2793        """
2794
2795        # Get interpolated values
2796        xmomentum = self.domain.get_quantity('xmomentum')
2797        ymomentum = self.domain.get_quantity('ymomentum')
2798
2799        uh = xmomentum.get_values(interpolation_points=self.midpoints,
2800                                  use_cache=True)
2801        vh = ymomentum.get_values(interpolation_points=self.midpoints,
2802                                  use_cache=True)
2803
2804        # Compute and sum flows across each segment
2805        total_flow = 0
2806        for i in range(len(uh)):
2807            # Inner product of momentum vector with segment normal [m^2/s]
2808            normal = self.segments[i].normal
2809            normal_momentum = uh[i]*normal[0] + vh[i]*normal[1]
2810
2811            # Flow across this segment [m^3/s]
2812            segment_flow = normal_momentum*self.segments[i].length
2813
2814            # Accumulate
2815            total_flow += segment_flow
2816
2817        return total_flow
2818 
2819
2820    ##
2821    # @brief calculate current energy flow through cross section
2822    def get_energy_through_cross_section(self):
2823        """Obtain average energy head [m] across specified cross section.
2824
2825        Output:
2826            E: Average energy [m] across given segments for all stored times.
2827
2828        The average velocity is computed for each triangle intersected by
2829        the polyline and averaged weighted by segment lengths.
2830
2831        The typical usage of this function would be to get average energy of
2832        flow in a channel, and the polyline would then be a cross section
2833        perpendicular to the flow.
2834
2835        #FIXME (Ole) - need name for this energy reflecting that its dimension
2836        is [m].
2837        """
2838
2839        # Get interpolated values
2840        stage = self.domain.get_quantity('stage')
2841        elevation = self.domain.get_quantity('elevation')
2842        xmomentum = self.domain.get_quantity('xmomentum')
2843        ymomentum = self.domain.get_quantity('ymomentum')
2844
2845        w = stage.get_values(interpolation_points=midpoints, use_cache=True)
2846        z = elevation.get_values(interpolation_points=midpoints, use_cache=True)
2847        uh = xmomentum.get_values(interpolation_points=midpoints,
2848                                  use_cache=True)
2849        vh = ymomentum.get_values(interpolation_points=midpoints,
2850                                  use_cache=True)
2851        h = w-z                # Depth
2852
2853        # Compute total length of polyline for use with weighted averages
2854        total_line_length = 0.0
2855        for segment in segments:
2856            total_line_length += segment.length
2857
2858        # Compute and sum flows across each segment
2859        average_energy = 0.0
2860        for i in range(len(w)):
2861            # Average velocity across this segment
2862            if h[i] > epsilon:
2863                # Use protection against degenerate velocities
2864                u = uh[i]/(h[i] + h0/h[i])
2865                v = vh[i]/(h[i] + h0/h[i])
2866            else:
2867                u = v = 0.0
2868
2869            speed_squared = u*u + v*v
2870            kinetic_energy = 0.5*speed_squared/g
2871
2872            if kind == 'specific':
2873                segment_energy = h[i] + kinetic_energy
2874            elif kind == 'total':
2875                segment_energy = w[i] + kinetic_energy
2876            else:
2877                msg = 'Energy kind must be either "specific" or "total".'
2878                msg += ' I got %s' %kind
2879
2880            # Add to weighted average
2881            weigth = segments[i].length/total_line_length
2882            average_energy += segment_energy*weigth
2883
2884        return average_energy
2885
2886
2887
2888################################################################################
2889# Initialise module
2890################################################################################
2891
2892from anuga.utilities import compile
2893if compile.can_use_C_extension('shallow_water_ext.c'):
2894    # Underlying C implementations can be accessed
2895    from shallow_water_ext import rotate, assign_windfield_values
2896else:
2897    msg = 'C implementations could not be accessed by %s.\n ' % __file__
2898    msg += 'Make sure compile_all.py has been run as described in '
2899    msg += 'the ANUGA installation guide.'
2900    raise Exception, msg
2901
2902# Optimisation with psyco
2903from anuga.config import use_psyco
2904if use_psyco:
2905    try:
2906        import psyco
2907    except:
2908        import os
2909        if os.name == 'posix' and os.uname()[4] in ['x86_64', 'ia64']:
2910            pass
2911            #Psyco isn't supported on 64 bit systems, but it doesn't matter
2912        else:
2913            msg = ('WARNING: psyco (speedup) could not be imported, '
2914                   'you may want to consider installing it')
2915            log.critical(msg)
2916    else:
2917        psyco.bind(Domain.distribute_to_vertices_and_edges)
2918        psyco.bind(Domain.compute_fluxes)
2919
2920
2921if __name__ == "__main__":
2922    pass
Note: See TracBrowser for help on using the repository browser.