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

Last change on this file since 7967 was 7967, checked in by steve, 13 years ago

Commits from session with Nariman

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