source: trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py @ 7810

Last change on this file since 7810 was 7810, checked in by hudson, 14 years ago

Added aggressive psyco optimisation, fixed benchmark app to work with new API.

  • Property svn:keywords set to LastChangedDate LastChangedRevision LastChangedBy HeadURL Id
File size: 49.7 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: trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py $
68
69Constraints: See GPL license in the user guide
70Version: 1.0 ($Revision: 7810 $)
71ModifiedBy:
72    $Author: hudson $
73    $Date: 2010-06-08 04:31:04 +0000 (Tue, 08 Jun 2010) $
74"""
75
76# Subversion keywords:
77#
78# $LastChangedDate: 2010-06-08 04:31:04 +0000 (Tue, 08 Jun 2010) $
79# $LastChangedRevision: 7810 $
80# $LastChangedBy: hudson $
81
82
83import numpy as num
84
85from anuga.abstract_2d_finite_volumes.generic_domain \
86                    import Generic_Domain
87
88from anuga.shallow_water.forcing import Cross_section
89from anuga.pmesh.mesh_interface import create_mesh_from_regions
90from anuga.utilities.numerical_tools import gradient, mean, ensure_numeric
91from anuga.geospatial_data.geospatial_data import ensure_geospatial
92from anuga.file.sww import SWW_file 
93from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
94
95from anuga.fit_interpolate.interpolate import Modeltime_too_late, \
96                                              Modeltime_too_early
97
98from anuga.geometry.polygon import inside_polygon, polygon_area, \
99                                    is_inside_polygon
100import anuga.utilities.log as log
101
102import types
103from types import IntType, FloatType
104from warnings import warn
105
106
107################################################################################
108# Shallow water domain
109################################################################################
110
111##
112# @brief Class for a shallow water domain.
113class Domain(Generic_Domain):
114
115    #conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
116    #other_quantities = ['elevation', 'friction']
117
118    ##
119    # @brief Instantiate a shallow water domain.
120    # @param coordinates
121    # @param vertices
122    # @param boundary
123    # @param tagged_elements
124    # @param geo_reference
125    # @param use_inscribed_circle
126    # @param mesh_filename
127    # @param use_cache
128    # @param verbose
129    # @param evolved_quantities
130    # @param full_send_dict
131    # @param ghost_recv_dict
132    # @param processor
133    # @param numproc
134    # @param number_of_full_nodes
135    # @param number_of_full_triangles
136    def __init__(self,
137                 coordinates=None,
138                 vertices=None,
139                 boundary=None,
140                 tagged_elements=None,
141                 geo_reference=None,
142                 use_inscribed_circle=False,
143                 mesh_filename=None,
144                 use_cache=False,
145                 verbose=False,
146                 conserved_quantities = None,
147                 evolved_quantities = None,
148                 other_quantities = None,
149                 full_send_dict=None,
150                 ghost_recv_dict=None,
151                 processor=0,
152                 numproc=1,
153                 number_of_full_nodes=None,
154                 number_of_full_triangles=None):
155
156        # Define quantities for the shallow_water domain
157        if conserved_quantities == None:
158            conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
159
160        if evolved_quantities == None:
161            evolved_quantities =  ['stage', 'xmomentum', 'ymomentum']
162           
163        if other_quantities == None:
164            other_quantities = ['elevation', 'friction']
165       
166        Generic_Domain.__init__(self,
167                                coordinates,
168                                vertices,
169                                boundary,
170                                conserved_quantities,
171                                evolved_quantities,
172                                other_quantities,
173                                tagged_elements,
174                                geo_reference,
175                                use_inscribed_circle,
176                                mesh_filename,
177                                use_cache,
178                                verbose,
179                                full_send_dict,
180                                ghost_recv_dict,
181                                processor,
182                                numproc,
183                                number_of_full_nodes=number_of_full_nodes,
184                                number_of_full_triangles=number_of_full_triangles)
185
186        self.set_defaults()
187
188 
189        self.forcing_terms.append(manning_friction_implicit)
190        self.forcing_terms.append(gravity)
191
192        # Stored output
193        self.store = True
194        self.set_store_vertices_uniquely(False)
195
196        self.quantities_to_be_stored = {'elevation': 1, 
197                                        'stage': 2, 
198                                        'xmomentum': 2, 
199                                        'ymomentum': 2}
200
201
202
203
204    ##
205    # @brief Set default values, usually read in from a config file
206    # @param flag
207    def set_defaults(self):
208        """Set the default values in this routine. That way we can inherit class
209        and just over redefine the defaults for the new class
210        """
211
212        from anuga.config import minimum_storable_height
213        from anuga.config import minimum_allowed_height, maximum_allowed_speed
214        from anuga.config import g, epsilon, beta_w, beta_w_dry,\
215             beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, tight_slope_limiters
216        from anuga.config import alpha_balance
217        from anuga.config import optimise_dry_cells
218        from anuga.config import optimised_gradient_limiter
219        from anuga.config import use_edge_limiter
220        from anuga.config import use_centroid_velocities
221
222
223
224        self.set_minimum_allowed_height(minimum_allowed_height)
225        self.maximum_allowed_speed = maximum_allowed_speed
226
227        self.g = g
228        self.beta_w = beta_w
229        self.beta_w_dry = beta_w_dry
230        self.beta_uh = beta_uh
231        self.beta_uh_dry = beta_uh_dry
232        self.beta_vh = beta_vh
233        self.beta_vh_dry = beta_vh_dry
234        self.alpha_balance = alpha_balance
235
236        self.tight_slope_limiters = tight_slope_limiters
237        self.optimise_dry_cells = optimise_dry_cells
238
239       
240        self.set_new_mannings_function(False)
241
242        self.minimum_storable_height = minimum_storable_height
243
244         # Limiters
245        self.use_edge_limiter = use_edge_limiter
246        self.optimised_gradient_limiter = optimised_gradient_limiter
247        self.use_centroid_velocities = use_centroid_velocities       
248
249    ##
250    # @brief
251    # @param flag
252    def set_new_mannings_function(self, flag=True):
253        """Cludge to allow unit test to pass, but to
254        also introduce new mannings friction function
255        which takes into account the slope of the bed.
256        The flag is tested in the python wrapper
257        mannings_friction_implicit
258        """
259        if flag:
260            self.use_new_mannings = True
261        else:
262            self.use_new_mannings = False
263
264
265    ##
266    # @brief
267    # @param flag
268    def set_use_edge_limiter(self, flag=True):
269        """Cludge to allow unit test to pass, but to
270        also introduce new edge limiting. The flag is
271        tested in distribute_to_vertices_and_edges
272        """
273        if flag:
274            self.use_edge_limiter = True
275        else:
276            self.use_edge_limiter = False
277
278
279         
280    ##
281    # @brief
282    # @param beta
283    def set_all_limiters(self, beta):
284        """Shorthand to assign one constant value [0,1] to all limiters.
285        0 Corresponds to first order, where as larger values make use of
286        the second order scheme.
287        """
288
289        self.beta_w = beta
290        self.beta_w_dry = beta
291        self.quantities['stage'].beta = beta
292
293        self.beta_uh = beta
294        self.beta_uh_dry = beta
295        self.quantities['xmomentum'].beta = beta
296
297        self.beta_vh = beta
298        self.beta_vh_dry = beta
299        self.quantities['ymomentum'].beta = beta
300
301    ##
302    # @brief
303    # @param flag
304    # @param reduction
305    def set_store_vertices_uniquely(self, flag, reduction=None):
306        """Decide whether vertex values should be stored uniquely as
307        computed in the model (True) or whether they should be reduced to one
308        value per vertex using self.reduction (False).
309        """
310
311        # FIXME (Ole): how about using the word "continuous vertex values" or
312        # "continuous stage surface"
313        self.smooth = not flag
314
315        # Reduction operation for get_vertex_values
316        if reduction is None:
317            self.reduction = mean
318            #self.reduction = min  #Looks better near steep slopes
319
320    ##
321    # @brief Set the minimum depth that will be written to an SWW file.
322    # @param minimum_storable_height The minimum stored height (in m).
323    def set_minimum_storable_height(self, minimum_storable_height):
324        """Set the minimum depth that will be recognised when writing
325        to an sww file. This is useful for removing thin water layers
326        that seems to be caused by friction creep.
327
328        The minimum allowed sww depth is in meters.
329        """
330
331        self.minimum_storable_height = minimum_storable_height
332
333    ##
334    # @brief
335    # @param minimum_allowed_height
336    def set_minimum_allowed_height(self, minimum_allowed_height):
337        """Set minimum depth that will be recognised in the numerical scheme.
338
339        The minimum allowed depth is in meters.
340
341        The parameter H0 (Minimal height for flux computation)
342        is also set by this function
343        """
344
345        #FIXME (Ole): rename H0 to minimum_allowed_height_in_flux_computation
346
347        #FIXME (Ole): Maybe use histogram to identify isolated extreme speeds
348        #and deal with them adaptively similarly to how we used to use 1 order
349        #steps to recover.
350
351        self.minimum_allowed_height = minimum_allowed_height
352        self.H0 = minimum_allowed_height
353
354    ##
355    # @brief
356    # @param maximum_allowed_speed
357    def set_maximum_allowed_speed(self, maximum_allowed_speed):
358        """Set the maximum particle speed that is allowed in water
359        shallower than minimum_allowed_height. This is useful for
360        controlling speeds in very thin layers of water and at the same time
361        allow some movement avoiding pooling of water.
362        """
363
364        self.maximum_allowed_speed = maximum_allowed_speed
365
366    ##
367    # @brief
368    # @param points_file_block_line_size
369    def set_points_file_block_line_size(self, points_file_block_line_size):
370        """Set the minimum depth that will be recognised when writing
371        to an sww file. This is useful for removing thin water layers
372        that seems to be caused by friction creep.
373
374        The minimum allowed sww depth is in meters.
375        """
376        self.points_file_block_line_size = points_file_block_line_size
377
378       
379    # FIXME: Probably obsolete in its curren form   
380    ##
381    # @brief Set the quantities that will be written to an SWW file.
382    # @param q The quantities to be written.
383    # @note Param 'q' may be None, single quantity or list of quantity strings.
384    # @note If 'q' is None, no quantities will be stored in the SWW file.
385    def set_quantities_to_be_stored(self, q):
386        """Specify which quantities will be stored in the sww file
387       
388        q must be either:
389          - a dictionary with quantity names
390          - a list of quantity names (for backwards compatibility)
391          - None
392
393        The format of the dictionary is as follows
394       
395        quantity_name: flag where flag must be either 1 or 2.
396        If flag is 1, the quantity is considered static and will be
397        stored once at the beginning of the simulation in a 1D array.
398       
399        If flag is 2, the quantity is considered time dependent and
400        it will be stored at each yieldstep by appending it to the
401        appropriate 2D array in the sww file.   
402       
403        If q is None, storage will be switched off altogether.
404       
405        Once the simulation has started and thw sww file opened,
406        this function will have no effect.
407       
408        The format, where q is a list of names is for backwards compatibility only.
409        It will take the specified quantities to be time dependent and assume
410        'elevation' to be static regardless.
411        """
412
413        if q is None:
414            self.quantities_to_be_stored = {}
415            self.store = False
416            return
417
418        # Check correcness
419        for quantity_name in q:
420            msg = ('Quantity %s is not a valid conserved quantity'
421                   % quantity_name)
422            assert quantity_name in self.quantities, msg
423
424        if type(q) == types.ListType:
425
426            msg = 'List arguments to set_quantities_to_be_stored '
427            msg += 'has been deprecated and will be removed in future '
428            msg += 'versions of ANUGA.'
429            msg += 'Please use dictionary argument instead'
430            warn(msg, DeprecationWarning) 
431
432       
433       
434            # FIXME: Raise deprecation warning
435            tmp = {}
436            for x in q:
437                tmp[x] = 2
438            tmp['elevation'] = 1   
439            q = tmp     
440           
441        assert type(q) == types.DictType   
442        self.quantities_to_be_stored = q
443
444    ##
445    # @brief
446    # @param indices
447    def get_wet_elements(self, indices=None):
448        """Return indices for elements where h > minimum_allowed_height
449
450        Optional argument:
451            indices is the set of element ids that the operation applies to.
452
453        Usage:
454            indices = get_wet_elements()
455
456        Note, centroid values are used for this operation
457        """
458
459        # Water depth below which it is considered to be 0 in the model
460        # FIXME (Ole): Allow this to be specified as a keyword argument as well
461        from anuga.config import minimum_allowed_height
462
463        elevation = self.get_quantity('elevation').\
464                        get_values(location='centroids', indices=indices)
465        stage = self.get_quantity('stage').\
466                    get_values(location='centroids', indices=indices)
467        depth = stage - elevation
468
469        # Select indices for which depth > 0
470        wet_indices = num.compress(depth > minimum_allowed_height,
471                                   num.arange(len(depth)))
472        return wet_indices
473
474    ##
475    # @brief
476    # @param indices
477    def get_maximum_inundation_elevation(self, indices=None):
478        """Return highest elevation where h > 0
479
480        Optional argument:
481            indices is the set of element ids that the operation applies to.
482
483        Usage:
484            q = get_maximum_inundation_elevation()
485
486        Note, centroid values are used for this operation
487        """
488
489        wet_elements = self.get_wet_elements(indices)
490        return self.get_quantity('elevation').\
491                   get_maximum_value(indices=wet_elements)
492
493    ##
494    # @brief
495    # @param indices
496    def get_maximum_inundation_location(self, indices=None):
497        """Return location of highest elevation where h > 0
498
499        Optional argument:
500            indices is the set of element ids that the operation applies to.
501
502        Usage:
503            q = get_maximum_inundation_location()
504
505        Note, centroid values are used for this operation
506        """
507
508        wet_elements = self.get_wet_elements(indices)
509        return self.get_quantity('elevation').\
510                   get_maximum_location(indices=wet_elements)
511
512
513    def get_flow_through_cross_section(self, polyline, verbose=False):
514        """Get the total flow through an arbitrary poly line.
515
516        This is a run-time equivalent of the function with same name
517        in sww_interrogate.py
518
519        Input:
520            polyline: Representation of desired cross section - it may contain
521                      multiple sections allowing for complex shapes. Assume
522                      absolute UTM coordinates.
523                      Format [[x0, y0], [x1, y1], ...]
524
525        Output:
526            Q: Total flow [m^3/s] across given segments.
527        """
528
529
530        cross_section = Cross_section(self, polyline, verbose)
531
532        return cross_section.get_flow_through_cross_section()
533
534
535    def get_energy_through_cross_section(self, polyline,
536                                         kind='total',
537                                         verbose=False):
538        """Obtain average energy head [m] across specified cross section.
539
540        Inputs:
541            polyline: Representation of desired cross section - it may contain
542                      multiple sections allowing for complex shapes. Assume
543                      absolute UTM coordinates.
544                      Format [[x0, y0], [x1, y1], ...]
545            kind:     Select which energy to compute.
546                      Options are 'specific' and 'total' (default)
547
548        Output:
549            E: Average energy [m] across given segments for all stored times.
550
551        The average velocity is computed for each triangle intersected by
552        the polyline and averaged weighted by segment lengths.
553
554        The typical usage of this function would be to get average energy of
555        flow in a channel, and the polyline would then be a cross section
556        perpendicular to the flow.
557
558        #FIXME (Ole) - need name for this energy reflecting that its dimension
559        is [m].
560        """
561
562
563
564        cross_section = Cross_section(self, polyline, verbose)
565
566        return cross_section.get_energy_through_cross_section(kind)
567
568
569    ##
570    # @brief Run integrity checks on shallow water domain.
571    def check_integrity(self):
572        Generic_Domain.check_integrity(self)
573
574        #Check that we are solving the shallow water wave equation
575        msg = 'First conserved quantity must be "stage"'
576        assert self.conserved_quantities[0] == 'stage', msg
577        msg = 'Second conserved quantity must be "xmomentum"'
578        assert self.conserved_quantities[1] == 'xmomentum', msg
579        msg = 'Third conserved quantity must be "ymomentum"'
580        assert self.conserved_quantities[2] == 'ymomentum', msg
581
582    ##
583    # @brief
584    def extrapolate_second_order_sw(self):
585        #Call correct module function (either from this module or C-extension)
586        extrapolate_second_order_sw(self)
587
588    ##
589    # @brief
590    def compute_fluxes(self):
591        #Call correct module function (either from this module or C-extension)
592        compute_fluxes(self)
593
594    ##
595    # @brief
596    def distribute_to_vertices_and_edges(self):
597        # Call correct module function
598        if self.use_edge_limiter:
599            distribute_using_edge_limiter(self)
600        else:
601            distribute_using_vertex_limiter(self)
602
603
604
605    ##
606    # @brief Evolve the model by one step.
607    # @param yieldstep
608    # @param finaltime
609    # @param duration
610    # @param skip_initial_step
611    def evolve(self,
612               yieldstep=None,
613               finaltime=None,
614               duration=None,
615               skip_initial_step=False):
616        """Specialisation of basic evolve method from parent class"""
617
618        # Call check integrity here rather than from user scripts
619        # self.check_integrity()
620
621        msg = 'Attribute self.beta_w must be in the interval [0, 2]'
622        assert 0 <= self.beta_w <= 2.0, msg
623
624        # Initial update of vertex and edge values before any STORAGE
625        # and or visualisation.
626        # This is done again in the initialisation of the Generic_Domain
627        # evolve loop but we do it here to ensure the values are ok for storage.
628        self.distribute_to_vertices_and_edges()
629
630        if self.store is True and self.time == 0.0:
631            self.initialise_storage()
632
633        # Call basic machinery from parent class
634        for t in Generic_Domain.evolve(self, yieldstep=yieldstep,
635                                       finaltime=finaltime, duration=duration,
636                                       skip_initial_step=skip_initial_step):
637            # Store model data, e.g. for subsequent visualisation
638            if self.store is True:
639                self.store_timestep()
640
641            # Pass control on to outer loop for more specific actions
642            yield(t)
643
644    ##
645    # @brief
646    def initialise_storage(self):
647        """Create and initialise self.writer object for storing data.
648        Also, save x,y and bed elevation
649        """
650       
651        # Initialise writer
652        self.writer = SWW_file(self)
653
654        # Store vertices and connectivity
655        self.writer.store_connectivity()
656
657    ##
658    # @brief
659    # @param name
660    def store_timestep(self):
661        """Store time dependent quantities and time.
662
663        Precondition:
664           self.writer has been initialised
665        """
666
667        self.writer.store_timestep()
668
669    ##
670    # @brief Get time stepping statistics string for printing.
671    # @param track_speeds
672    # @param triangle_id
673    def timestepping_statistics(self,
674                                track_speeds=False,
675                                triangle_id=None):
676        """Return string with time stepping statistics for printing or logging
677
678        Optional boolean keyword track_speeds decides whether to report
679        location of smallest timestep as well as a histogram and percentile
680        report.
681        """
682
683        from anuga.config import epsilon, g
684
685        # Call basic machinery from parent class
686        msg = Generic_Domain.timestepping_statistics(self, track_speeds,
687                                                     triangle_id)
688
689        if track_speeds is True:
690            # qwidth determines the text field used for quantities
691            qwidth = self.qwidth
692
693            # Selected triangle
694            k = self.k
695
696            # Report some derived quantities at vertices, edges and centroid
697            # specific to the shallow water wave equation
698            z = self.quantities['elevation']
699            w = self.quantities['stage']
700
701            Vw = w.get_values(location='vertices', indices=[k])[0]
702            Ew = w.get_values(location='edges', indices=[k])[0]
703            Cw = w.get_values(location='centroids', indices=[k])
704
705            Vz = z.get_values(location='vertices', indices=[k])[0]
706            Ez = z.get_values(location='edges', indices=[k])[0]
707            Cz = z.get_values(location='centroids', indices=[k])
708
709            name = 'depth'
710            Vh = Vw-Vz
711            Eh = Ew-Ez
712            Ch = Cw-Cz
713
714            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
715                 %(name.ljust(qwidth), Vh[0], Vh[1], Vh[2])
716
717            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
718                 %(name.ljust(qwidth), Eh[0], Eh[1], Eh[2])
719
720            s += '    %s: centroid_value = %.4f\n'\
721                 %(name.ljust(qwidth), Ch[0])
722
723            msg += s
724
725            uh = self.quantities['xmomentum']
726            vh = self.quantities['ymomentum']
727
728            Vuh = uh.get_values(location='vertices', indices=[k])[0]
729            Euh = uh.get_values(location='edges', indices=[k])[0]
730            Cuh = uh.get_values(location='centroids', indices=[k])
731
732            Vvh = vh.get_values(location='vertices', indices=[k])[0]
733            Evh = vh.get_values(location='edges', indices=[k])[0]
734            Cvh = vh.get_values(location='centroids', indices=[k])
735
736            # Speeds in each direction
737            Vu = Vuh/(Vh + epsilon)
738            Eu = Euh/(Eh + epsilon)
739            Cu = Cuh/(Ch + epsilon)
740            name = 'U'
741            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
742                 %(name.ljust(qwidth), Vu[0], Vu[1], Vu[2])
743
744            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
745                 %(name.ljust(qwidth), Eu[0], Eu[1], Eu[2])
746
747            s += '    %s: centroid_value = %.4f\n'\
748                 %(name.ljust(qwidth), Cu[0])
749
750            msg += s
751
752            Vv = Vvh/(Vh + epsilon)
753            Ev = Evh/(Eh + epsilon)
754            Cv = Cvh/(Ch + epsilon)
755            name = 'V'
756            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
757                 %(name.ljust(qwidth), Vv[0], Vv[1], Vv[2])
758
759            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
760                 %(name.ljust(qwidth), Ev[0], Ev[1], Ev[2])
761
762            s += '    %s: centroid_value = %.4f\n'\
763                 %(name.ljust(qwidth), Cv[0])
764
765            msg += s
766
767            # Froude number in each direction
768            name = 'Froude (x)'
769            Vfx = Vu/(num.sqrt(g*Vh) + epsilon)
770            Efx = Eu/(num.sqrt(g*Eh) + epsilon)
771            Cfx = Cu/(num.sqrt(g*Ch) + epsilon)
772
773            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
774                 %(name.ljust(qwidth), Vfx[0], Vfx[1], Vfx[2])
775
776            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
777                 %(name.ljust(qwidth), Efx[0], Efx[1], Efx[2])
778
779            s += '    %s: centroid_value = %.4f\n'\
780                 %(name.ljust(qwidth), Cfx[0])
781
782            msg += s
783
784            name = 'Froude (y)'
785            Vfy = Vv/(num.sqrt(g*Vh) + epsilon)
786            Efy = Ev/(num.sqrt(g*Eh) + epsilon)
787            Cfy = Cv/(num.sqrt(g*Ch) + epsilon)
788
789            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
790                 %(name.ljust(qwidth), Vfy[0], Vfy[1], Vfy[2])
791
792            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
793                 %(name.ljust(qwidth), Efy[0], Efy[1], Efy[2])
794
795            s += '    %s: centroid_value = %.4f\n'\
796                 %(name.ljust(qwidth), Cfy[0])
797
798            msg += s
799
800        return msg
801       
802       
803
804    def compute_boundary_flows(self):
805        """Compute boundary flows at current timestep.
806                       
807        Quantities computed are:
808           Total inflow across boundary
809           Total outflow across boundary
810           Flow across each tagged boundary segment
811        """
812               
813        # Run through boundary array and compute for each segment
814        # the normal momentum ((uh, vh) dot normal) times segment length.
815        # Based on sign accumulate this into boundary_inflow and boundary_outflow.
816                       
817        # Compute flows along boundary
818       
819        uh = self.get_quantity('xmomentum').get_values(location='edges')
820        vh = self.get_quantity('ymomentum').get_values(location='edges')       
821       
822        # Loop through edges that lie on the boundary and calculate
823        # flows
824        boundary_flows = {}
825        total_boundary_inflow = 0.0
826        total_boundary_outflow = 0.0
827        for vol_id, edge_id in self.boundary:
828            # Compute normal flow across edge. Since normal vector points
829            # away from triangle, a positive sign means that water
830            # flows *out* from this triangle.
831             
832            momentum = [uh[vol_id, edge_id], vh[vol_id, edge_id]]
833            normal = self.mesh.get_normal(vol_id, edge_id)
834            length = self.mesh.get_edgelength(vol_id, edge_id)           
835            normal_flow = num.dot(momentum, normal)*length
836           
837            # Reverse sign so that + is taken to mean inflow
838            # and - means outflow. This is more intuitive.
839            edge_flow = -normal_flow
840           
841            # Tally up inflows and outflows separately
842            if edge_flow > 0:
843                # Flow is inflow     
844                total_boundary_inflow += edge_flow                                 
845            else:
846                # Flow is outflow
847                total_boundary_outflow += edge_flow   
848
849            # Tally up flows by boundary tag
850            tag = self.boundary[(vol_id, edge_id)]
851           
852            if tag not in boundary_flows:
853                boundary_flows[tag] = 0.0
854            boundary_flows[tag] += edge_flow
855           
856               
857        return boundary_flows, total_boundary_inflow, total_boundary_outflow
858       
859
860    def compute_forcing_flows(self):
861        """Compute flows in and out of domain due to forcing terms.
862                       
863        Quantities computed are:
864               
865       
866           Total inflow through forcing terms
867           Total outflow through forcing terms
868           Current total volume in domain       
869
870        """
871
872        #FIXME(Ole): We need to separate what part of explicit_update was
873        # due to the normal flux calculations and what is due to forcing terms.
874       
875        pass
876                       
877       
878    def compute_total_volume(self):
879        """Compute total volume (m^3) of water in entire domain
880        """
881       
882        area = self.mesh.get_areas()
883        volume = 0.0
884       
885        stage = self.get_quantity('stage').get_values(location='centroids')
886        elevation = self.get_quantity('elevation').get_values(location='centroids')       
887        depth = stage-elevation
888       
889        return num.sum(depth*area)
890       
891       
892    def volumetric_balance_statistics(self):               
893        """Create volumetric balance report suitable for printing or logging.
894        """
895       
896        (boundary_flows, total_boundary_inflow,
897         total_boundary_outflow) = self.compute_boundary_flows() 
898       
899        s = '---------------------------\n'       
900        s += 'Volumetric balance report:\n'
901        s += '--------------------------\n'
902        s += 'Total boundary inflow [m^3/s]: %.2f\n' % total_boundary_inflow
903        s += 'Total boundary outflow [m^3/s]: %.2f\n' % total_boundary_outflow       
904        s += 'Net boundary flow by tags [m^3/s]\n'
905        for tag in boundary_flows:
906            s += '    %s [m^3/s]: %.2f\n' % (tag, boundary_flows[tag])
907       
908        s += 'Total net boundary flow [m^3/s]: %.2f\n' % (total_boundary_inflow + total_boundary_outflow) 
909        s += 'Total volume in domain [m^3]: %.2f\n' % self.compute_total_volume()
910       
911        # The go through explicit forcing update and record the rate of change for stage and
912        # record into forcing_inflow and forcing_outflow. Finally compute integral
913        # of depth to obtain total volume of domain.
914       
915        # FIXME(Ole): This part is not yet done.               
916       
917        return s       
918           
919################################################################################
920# End of class Shallow Water Domain
921################################################################################
922
923#-----------------
924# Flux computation
925#-----------------
926
927## @brief Compute fluxes and timestep suitable for all volumes in domain.
928# @param domain The domain to calculate fluxes for.
929def compute_fluxes(domain):
930    """Compute fluxes and timestep suitable for all volumes in domain.
931
932    Compute total flux for each conserved quantity using "flux_function"
933
934    Fluxes across each edge are scaled by edgelengths and summed up
935    Resulting flux is then scaled by area and stored in
936    explicit_update for each of the three conserved quantities
937    stage, xmomentum and ymomentum
938
939    The maximal allowable speed computed by the flux_function for each volume
940    is converted to a timestep that must not be exceeded. The minimum of
941    those is computed as the next overall timestep.
942
943    Post conditions:
944      domain.explicit_update is reset to computed flux values
945      domain.timestep is set to the largest step satisfying all volumes.
946
947    This wrapper calls the underlying C version of compute fluxes
948    """
949
950    import sys
951    from shallow_water_ext import compute_fluxes_ext_central \
952                                  as compute_fluxes_ext
953
954    N = len(domain)    # number_of_triangles
955
956    # Shortcuts
957    Stage = domain.quantities['stage']
958    Xmom = domain.quantities['xmomentum']
959    Ymom = domain.quantities['ymomentum']
960    Bed = domain.quantities['elevation']
961
962    timestep = float(sys.maxint)
963
964    flux_timestep = compute_fluxes_ext(timestep,
965                                       domain.epsilon,
966                                       domain.H0,
967                                       domain.g,
968                                       domain.neighbours,
969                                       domain.neighbour_edges,
970                                       domain.normals,
971                                       domain.edgelengths,
972                                       domain.radii,
973                                       domain.areas,
974                                       domain.tri_full_flag,
975                                       Stage.edge_values,
976                                       Xmom.edge_values,
977                                       Ymom.edge_values,
978                                       Bed.edge_values,
979                                       Stage.boundary_values,
980                                       Xmom.boundary_values,
981                                       Ymom.boundary_values,
982                                       Stage.explicit_update,
983                                       Xmom.explicit_update,
984                                       Ymom.explicit_update,
985                                       domain.already_computed_flux,
986                                       domain.max_speed,
987                                       int(domain.optimise_dry_cells))
988
989    domain.flux_timestep = flux_timestep
990
991################################################################################
992# Module functions for gradient limiting
993################################################################################
994
995##
996# @brief Wrapper for C version of extrapolate_second_order_sw.
997# @param domain The domain to operate on.
998# @note MH090605 The following method belongs to the shallow_water domain class
999#       see comments in the corresponding method in shallow_water_ext.c
1000def extrapolate_second_order_sw(domain):
1001    """Wrapper calling C version of extrapolate_second_order_sw"""
1002
1003    import sys
1004    from shallow_water_ext import extrapolate_second_order_sw as extrapol2
1005
1006    N = len(domain) # number_of_triangles
1007
1008    # Shortcuts
1009    Stage = domain.quantities['stage']
1010    Xmom = domain.quantities['xmomentum']
1011    Ymom = domain.quantities['ymomentum']
1012    Elevation = domain.quantities['elevation']
1013
1014    extrapol2(domain,
1015              domain.surrogate_neighbours,
1016              domain.number_of_boundaries,
1017              domain.centroid_coordinates,
1018              Stage.centroid_values,
1019              Xmom.centroid_values,
1020              Ymom.centroid_values,
1021              Elevation.centroid_values,
1022              domain.vertex_coordinates,
1023              Stage.vertex_values,
1024              Xmom.vertex_values,
1025              Ymom.vertex_values,
1026              Elevation.vertex_values,
1027              int(domain.optimise_dry_cells))
1028
1029##
1030# @brief Distribution from centroids to vertices specific to the SWW eqn.
1031# @param domain The domain to operate on.
1032def distribute_using_vertex_limiter(domain):
1033    """Distribution from centroids to vertices specific to the SWW equation.
1034
1035    It will ensure that h (w-z) is always non-negative even in the
1036    presence of steep bed-slopes by taking a weighted average between shallow
1037    and deep cases.
1038
1039    In addition, all conserved quantities get distributed as per either a
1040    constant (order==1) or a piecewise linear function (order==2).
1041
1042    FIXME: more explanation about removal of artificial variability etc
1043
1044    Precondition:
1045      All quantities defined at centroids and bed elevation defined at
1046      vertices.
1047
1048    Postcondition
1049      Conserved quantities defined at vertices
1050    """
1051
1052    # Remove very thin layers of water
1053    protect_against_infinitesimal_and_negative_heights(domain)
1054
1055    # Extrapolate all conserved quantities
1056    if domain.optimised_gradient_limiter:
1057        # MH090605 if second order,
1058        # perform the extrapolation and limiting on
1059        # all of the conserved quantities
1060
1061        if (domain._order_ == 1):
1062            for name in domain.conserved_quantities:
1063                Q = domain.quantities[name]
1064                Q.extrapolate_first_order()
1065        elif domain._order_ == 2:
1066            domain.extrapolate_second_order_sw()
1067        else:
1068            raise 'Unknown order'
1069    else:
1070        # Old code:
1071        for name in domain.conserved_quantities:
1072            Q = domain.quantities[name]
1073
1074            if domain._order_ == 1:
1075                Q.extrapolate_first_order()
1076            elif domain._order_ == 2:
1077                Q.extrapolate_second_order_and_limit_by_vertex()
1078            else:
1079                raise 'Unknown order'
1080
1081    # Take bed elevation into account when water heights are small
1082    balance_deep_and_shallow(domain)
1083
1084    # Compute edge values by interpolation
1085    for name in domain.conserved_quantities:
1086        Q = domain.quantities[name]
1087        Q.interpolate_from_vertices_to_edges()
1088
1089##
1090# @brief Distribution from centroids to edges specific to the SWW eqn.
1091# @param domain The domain to operate on.
1092def distribute_using_edge_limiter(domain):
1093    """Distribution from centroids to edges specific to the SWW eqn.
1094
1095    It will ensure that h (w-z) is always non-negative even in the
1096    presence of steep bed-slopes by taking a weighted average between shallow
1097    and deep cases.
1098
1099    In addition, all conserved quantities get distributed as per either a
1100    constant (order==1) or a piecewise linear function (order==2).
1101
1102
1103    Precondition:
1104      All quantities defined at centroids and bed elevation defined at
1105      vertices.
1106
1107    Postcondition
1108      Conserved quantities defined at vertices
1109    """
1110
1111    # Remove very thin layers of water
1112    protect_against_infinitesimal_and_negative_heights(domain)
1113
1114    for name in domain.conserved_quantities:
1115        Q = domain.quantities[name]
1116        if domain._order_ == 1:
1117            Q.extrapolate_first_order()
1118        elif domain._order_ == 2:
1119            Q.extrapolate_second_order_and_limit_by_edge()
1120        else:
1121            raise 'Unknown order'
1122
1123    balance_deep_and_shallow(domain)
1124
1125    # Compute edge values by interpolation
1126    for name in domain.conserved_quantities:
1127        Q = domain.quantities[name]
1128        Q.interpolate_from_vertices_to_edges()
1129
1130##
1131# @brief  Protect against infinitesimal heights and associated high velocities.
1132# @param domain The domain to operate on.
1133def protect_against_infinitesimal_and_negative_heights(domain):
1134    """Protect against infinitesimal heights and associated high velocities"""
1135
1136    from shallow_water_ext import protect
1137
1138    # Shortcuts
1139    wc = domain.quantities['stage'].centroid_values
1140    zc = domain.quantities['elevation'].centroid_values
1141    xmomc = domain.quantities['xmomentum'].centroid_values
1142    ymomc = domain.quantities['ymomentum'].centroid_values
1143
1144    protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
1145            domain.epsilon, wc, zc, xmomc, ymomc)
1146
1147##
1148# @brief Wrapper for C function balance_deep_and_shallow_c().
1149# @param domain The domain to operate on.
1150def balance_deep_and_shallow(domain):
1151    """Compute linear combination between stage as computed by
1152    gradient-limiters limiting using w, and stage computed by
1153    gradient-limiters limiting using h (h-limiter).
1154    The former takes precedence when heights are large compared to the
1155    bed slope while the latter takes precedence when heights are
1156    relatively small.  Anything in between is computed as a balanced
1157    linear combination in order to avoid numerical disturbances which
1158    would otherwise appear as a result of hard switching between
1159    modes.
1160
1161    Wrapper for C implementation
1162    """
1163
1164    from shallow_water_ext import balance_deep_and_shallow \
1165                                  as balance_deep_and_shallow_c
1166
1167    # Shortcuts
1168    wc = domain.quantities['stage'].centroid_values
1169    zc = domain.quantities['elevation'].centroid_values
1170    wv = domain.quantities['stage'].vertex_values
1171    zv = domain.quantities['elevation'].vertex_values
1172
1173    # Momentums at centroids
1174    xmomc = domain.quantities['xmomentum'].centroid_values
1175    ymomc = domain.quantities['ymomentum'].centroid_values
1176
1177    # Momentums at vertices
1178    xmomv = domain.quantities['xmomentum'].vertex_values
1179    ymomv = domain.quantities['ymomentum'].vertex_values
1180
1181    balance_deep_and_shallow_c(domain,
1182                               wc, zc, wv, zv, wc,
1183                               xmomc, ymomc, xmomv, ymomv)
1184
1185
1186
1187################################################################################
1188# Standard forcing terms
1189################################################################################
1190
1191##
1192# @brief Apply gravitational pull in the presence of bed slope.
1193# @param domain The domain to apply gravity to.
1194# @note Wrapper for C function gravity_c().
1195def gravity(domain):
1196    """Apply gravitational pull in the presence of bed slope
1197    Wrapper calls underlying C implementation
1198    """
1199
1200    from shallow_water_ext import gravity as gravity_c
1201
1202    xmom_update = domain.quantities['xmomentum'].explicit_update
1203    ymom_update = domain.quantities['ymomentum'].explicit_update
1204
1205    stage = domain.quantities['stage']
1206    elevation = domain.quantities['elevation']
1207
1208    h = stage.centroid_values - elevation.centroid_values
1209    z = elevation.vertex_values
1210
1211    x = domain.get_vertex_coordinates()
1212    g = domain.g
1213
1214    gravity_c(g, h, z, x, xmom_update, ymom_update)    #, 1.0e-6)
1215
1216##
1217# @brief Apply friction to a surface (implicit).
1218# @param domain The domain to apply Manning friction to.
1219# @note Wrapper for C function manning_friction_c().
1220def manning_friction_implicit(domain):
1221    """Apply (Manning) friction to water momentum
1222    Wrapper for c version
1223    """
1224
1225    from shallow_water_ext import manning_friction_old
1226    from shallow_water_ext import manning_friction_new
1227
1228    xmom = domain.quantities['xmomentum']
1229    ymom = domain.quantities['ymomentum']
1230
1231    x = domain.get_vertex_coordinates()
1232   
1233    w = domain.quantities['stage'].centroid_values
1234    z = domain.quantities['elevation'].vertex_values
1235
1236    uh = xmom.centroid_values
1237    vh = ymom.centroid_values
1238    eta = domain.quantities['friction'].centroid_values
1239
1240    xmom_update = xmom.semi_implicit_update
1241    ymom_update = ymom.semi_implicit_update
1242
1243    N = len(domain)
1244    eps = domain.minimum_allowed_height
1245    g = domain.g
1246
1247    if domain.use_new_mannings:
1248        manning_friction_new(g, eps, x, w, uh, vh, z, eta, xmom_update, ymom_update)
1249    else:
1250        manning_friction_old(g, eps, w, uh, vh, z, eta, xmom_update, ymom_update)
1251   
1252
1253##
1254# @brief Apply friction to a surface (explicit).
1255# @param domain The domain to apply Manning friction to.
1256# @note Wrapper for C function manning_friction_c().
1257def manning_friction_explicit(domain):
1258    """Apply (Manning) friction to water momentum
1259    Wrapper for c version
1260    """
1261
1262    from shallow_water_ext import manning_friction_old
1263    from shallow_water_ext import manning_friction_new
1264
1265    xmom = domain.quantities['xmomentum']
1266    ymom = domain.quantities['ymomentum']
1267
1268    x = domain.get_vertex_coordinates()
1269   
1270    w = domain.quantities['stage'].centroid_values
1271    z = domain.quantities['elevation'].vertex_values
1272
1273    uh = xmom.centroid_values
1274    vh = ymom.centroid_values
1275    eta = domain.quantities['friction'].centroid_values
1276
1277    xmom_update = xmom.explicit_update
1278    ymom_update = ymom.explicit_update
1279
1280    N = len(domain)
1281    eps = domain.minimum_allowed_height
1282    g = domain.g
1283
1284
1285    if domain.use_new_mannings:
1286        manning_friction_new(g, eps, x, w, uh, vh, z, eta, xmom_update, ymom_update)
1287    else:
1288        manning_friction_old(g, eps, w, uh, vh, z, eta, xmom_update, ymom_update)
1289
1290
1291
1292# FIXME (Ole): This was implemented for use with one of the analytical solutions (Sampson?)
1293##
1294# @brief Apply linear friction to a surface.
1295# @param domain The domain to apply Manning friction to.
1296# @note Is this still used (30 Oct 2007)?
1297def linear_friction(domain):
1298    """Apply linear friction to water momentum
1299
1300    Assumes quantity: 'linear_friction' to be present
1301    """
1302
1303    from math import sqrt
1304
1305    w = domain.quantities['stage'].centroid_values
1306    z = domain.quantities['elevation'].centroid_values
1307    h = w-z
1308
1309    uh = domain.quantities['xmomentum'].centroid_values
1310    vh = domain.quantities['ymomentum'].centroid_values
1311    tau = domain.quantities['linear_friction'].centroid_values
1312
1313    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1314    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1315
1316    N = len(domain) # number_of_triangles
1317    eps = domain.minimum_allowed_height
1318    g = domain.g #Not necessary? Why was this added?
1319
1320    for k in range(N):
1321        if tau[k] >= eps:
1322            if h[k] >= eps:
1323                S = -tau[k]/h[k]
1324
1325                #Update momentum
1326                xmom_update[k] += S*uh[k]
1327                ymom_update[k] += S*vh[k]
1328
1329def depth_dependent_friction(domain, default_friction,
1330                             surface_roughness_data,
1331                             verbose=False):
1332    """Returns an array of friction values for each wet element adjusted for depth.
1333
1334    Inputs:
1335        domain - computational domain object
1336        default_friction - depth independent bottom friction
1337        surface_roughness_data - N x 5 array of n0, d1, n1, d2, n2 values for each
1338        friction region.
1339
1340    Outputs:
1341        wet_friction - Array that can be used directly to update friction as follows:
1342                       domain.set_quantity('friction', wet_friction)
1343
1344       
1345       
1346    """
1347
1348    import numpy as num
1349   
1350    # Create a temp array to store updated depth dependent friction for wet elements
1351    # EHR this is outwardly inneficient but not obvious how to avoid recreating each call??????
1352    N=len(domain)
1353    wet_friction    = num.zeros(N, num.float)
1354    wet_friction[:] = default_n0   # Initially assign default_n0 to all array so sure have no zeros values
1355   
1356   
1357    depth = domain.create_quantity_from_expression('stage - elevation')  # create depth instance for this timestep
1358    # Recompute depth as vector 
1359    d = depth.get_values(location='centroids')
1360 
1361    # rebuild the 'friction' values adjusted for depth at this instant
1362    for i in domain.get_wet_elements():                                  # loop for each wet element in domain
1363       
1364        # Get roughness data for each element
1365        n0 = float(surface_roughness_data[i,0])
1366        d1 = float(surface_roughness_data[i,1])
1367        n1 = float(surface_roughness_data[i,2])
1368        d2 = float(surface_roughness_data[i,3])
1369        n2 = float(surface_roughness_data[i,4])
1370       
1371       
1372        # Recompute friction values from depth for this element
1373               
1374        if d[i]   <= d1: 
1375            depth_dependent_friction = n1
1376        elif d[i] >= d2:
1377            depth_dependent_friction = n2
1378        else:
1379            depth_dependent_friction = n1+((n2-n1)/(d2-d1))*(d[i]-d1)
1380           
1381        # check sanity of result
1382        if (depth_dependent_friction  < 0.010 or depth_dependent_friction > 9999.0) :
1383            log.critical('%s >>>> WARNING: computed depth_dependent friction '
1384                         'out of range, ddf%f, n1=%f, n2=%f'
1385                         % (model_data.basename,
1386                            depth_dependent_friction, n1, n2))
1387       
1388        # update depth dependent friction  for that wet element
1389        wet_friction[i] = depth_dependent_friction
1390       
1391    # EHR add code to show range of 'friction across domain at this instant as sanity check?????????
1392   
1393    if verbose :
1394        nvals=domain.get_quantity('friction').get_values(location='centroids')        # return array of domain nvals
1395        n_min=min(nvals)
1396        n_max=max(nvals)
1397       
1398        log.critical('         ++++ calculate_depth_dependent_friction - '
1399                     'Updated friction - range  %7.3f to %7.3f'
1400                     % (n_min, n_max))
1401   
1402    return wet_friction
1403
1404
1405
1406################################################################################
1407# Initialise module
1408################################################################################
1409
1410from anuga.utilities import compile
1411if compile.can_use_C_extension('shallow_water_ext.c'):
1412    # Underlying C implementations can be accessed
1413    from shallow_water_ext import assign_windfield_values
1414else:
1415    msg = 'C implementations could not be accessed by %s.\n ' % __file__
1416    msg += 'Make sure compile_all.py has been run as described in '
1417    msg += 'the ANUGA installation guide.'
1418    raise Exception, msg
1419
1420# Optimisation with psyco
1421#from anuga.config import use_psyco
1422#if use_psyco:
1423    #try:
1424        #import psyco
1425    #except:
1426        #import os
1427        #if os.name == 'posix' and os.uname()[4] in ['x86_64', 'ia64']:
1428            #pass
1429            ##Psyco isn't supported on 64 bit systems, but it doesn't matter
1430        #else:
1431            #msg = ('WARNING: psyco (speedup) could not be imported, '
1432                   #'you may want to consider installing it')
1433            #log.critical(msg)
1434    #else:
1435        #psyco.bind(Domain.distribute_to_vertices_and_edges)
1436        #psyco.bind(Domain.compute_fluxes)
1437
1438
1439if __name__ == "__main__":
1440    pass
Note: See TracBrowser for help on using the repository browser.