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

Last change on this file since 8124 was 8124, checked in by wilsonr, 13 years ago

Made changes described in ticket 359.

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