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

Last change on this file since 7731 was 7731, checked in by hudson, 13 years ago

Split boundaries out of shallow_water_domain module.

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