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

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

Adding in culvert routines into structures directory

  • Property svn:keywords set to LastChangedDate LastChangedRevision LastChangedBy HeadURL Id
File size: 47.7 KB
Line 
1"""
2Finite-volume computations of the shallow water wave equation.
3
4Title: ANGUA shallow_water_domain - 2D triangular domains for finite-volume
5       computations of the shallow water wave equation.
6
7
8Author: Ole Nielsen, Ole.Nielsen@ga.gov.au
9        Stephen Roberts, Stephen.Roberts@anu.edu.au
10        Duncan Gray, Duncan.Gray@ga.gov.au
11
12CreationDate: 2004
13
14Description:
15    This module contains a specialisation of class Domain from
16    module domain.py consisting of methods specific to the
17    Shallow Water Wave Equation
18
19    U_t + E_x + G_y = S
20
21    where
22
23    U = [w, uh, vh]
24    E = [uh, u^2h + gh^2/2, uvh]
25    G = [vh, uvh, v^2h + gh^2/2]
26    S represents source terms forcing the system
27    (e.g. gravity, friction, wind stress, ...)
28
29    and _t, _x, _y denote the derivative with respect to t, x and y
30    respectively.
31
32
33    The quantities are
34
35    symbol    variable name    explanation
36    x         x                horizontal distance from origin [m]
37    y         y                vertical distance from origin [m]
38    z         elevation        elevation of bed on which flow is modelled [m]
39    h         height           water height above z [m]
40    w         stage            absolute water level, w = z+h [m]
41    u                          speed in the x direction [m/s]
42    v                          speed in the y direction [m/s]
43    uh        xmomentum        momentum in the x direction [m^2/s]
44    vh        ymomentum        momentum in the y direction [m^2/s]
45
46    eta                        mannings friction coefficient [to appear]
47    nu                         wind stress coefficient [to appear]
48
49    The conserved quantities are w, uh, vh
50
51Reference:
52    Catastrophic Collapse of Water Supply Reservoirs in Urban Areas,
53    Christopher Zoppou and Stephen Roberts,
54    Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999
55
56    Hydrodynamic modelling of coastal inundation.
57    Nielsen, O., S. Roberts, D. Gray, A. McPherson and A. Hitchman
58    In Zerger, A. and Argent, R.M. (eds) MODSIM 2005 International Congress on
59    Modelling and Simulation. Modelling and Simulation Society of Australia and
60    New Zealand, December 2005, pp. 518-523. ISBN: 0-9758400-2-9.
61    http://www.mssanz.org.au/modsim05/papers/nielsen.pdf
62
63
64SeeAlso:
65    TRAC administration of ANUGA (User Manuals etc) at
66    https://datamining.anu.edu.au/anuga and Subversion repository at
67    $HeadURL: 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: 7939 $)
72ModifiedBy:
73    $Author: steve $
74    $Date: 2010-08-13 05:35:39 +0000 (Fri, 13 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 apply_fractional_steps(self):
536        """Loop throughand apply all extra fractional steps"""
537        for operator in self.fractional_step_operators:
538            operator.apply()
539
540    def distribute_to_vertices_and_edges(self):
541        """ Call correct module function """
542        if self.use_edge_limiter:
543            distribute_using_edge_limiter(self)
544        else:
545            distribute_using_vertex_limiter(self)
546
547
548
549    def evolve(self,
550               yieldstep=None,
551               finaltime=None,
552               duration=None,
553               skip_initial_step=False):
554        """Specialisation of basic evolve method from parent class.
555       
556            Evolve the model by 1 step.
557        """
558
559        # Call check integrity here rather than from user scripts
560        # self.check_integrity()
561
562        msg = 'Attribute self.beta_w must be in the interval [0, 2]'
563        assert 0 <= self.beta_w <= 2.0, msg
564
565        # Initial update of vertex and edge values before any STORAGE
566        # and or visualisation.
567        # This is done again in the initialisation of the Generic_Domain
568        # evolve loop but we do it here to ensure the values are ok for storage.
569        self.distribute_to_vertices_and_edges()
570
571        if self.store is True and self.time == 0.0:
572            self.initialise_storage()
573
574        # Call basic machinery from parent class
575        for t in Generic_Domain.evolve(self, yieldstep=yieldstep,
576                                       finaltime=finaltime, duration=duration,
577                                       skip_initial_step=skip_initial_step):
578            # Store model data, e.g. for subsequent visualisation
579            if self.store is True:
580                self.store_timestep()
581
582            # Pass control on to outer loop for more specific actions
583            yield(t)
584
585
586    def initialise_storage(self):
587        """Create and initialise self.writer object for storing data.
588        Also, save x,y and bed elevation
589        """
590       
591        # Initialise writer
592        self.writer = SWW_file(self)
593
594        # Store vertices and connectivity
595        self.writer.store_connectivity()
596
597
598    def store_timestep(self):
599        """Store time dependent quantities and time.
600
601        Precondition:
602           self.writer has been initialised
603        """
604
605        self.writer.store_timestep()
606
607
608    def timestepping_statistics(self,
609                                track_speeds=False,
610                                triangle_id=None):
611        """Return string with time stepping statistics for printing or logging
612
613        Optional boolean keyword track_speeds decides whether to report
614        location of smallest timestep as well as a histogram and percentile
615        report.
616        """
617
618        from anuga.config import epsilon, g
619
620        # Call basic machinery from parent class
621        msg = Generic_Domain.timestepping_statistics(self, track_speeds,
622                                                     triangle_id)
623
624        if track_speeds is True:
625            # qwidth determines the text field used for quantities
626            qwidth = self.qwidth
627
628            # Selected triangle
629            k = self.k
630
631            # Report some derived quantities at vertices, edges and centroid
632            # specific to the shallow water wave equation
633            z = self.quantities['elevation']
634            w = self.quantities['stage']
635
636            Vw = w.get_values(location='vertices', indices=[k])[0]
637            Ew = w.get_values(location='edges', indices=[k])[0]
638            Cw = w.get_values(location='centroids', indices=[k])
639
640            Vz = z.get_values(location='vertices', indices=[k])[0]
641            Ez = z.get_values(location='edges', indices=[k])[0]
642            Cz = z.get_values(location='centroids', indices=[k])
643
644            name = 'depth'
645            Vh = Vw-Vz
646            Eh = Ew-Ez
647            Ch = Cw-Cz
648
649            message  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
650                 % (name.ljust(qwidth), Vh[0], Vh[1], Vh[2])
651
652            message += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
653                 % (name.ljust(qwidth), Eh[0], Eh[1], Eh[2])
654
655            message += '    %s: centroid_value = %.4f\n'\
656                 % (name.ljust(qwidth), Ch[0])
657
658            msg += message
659
660            uh = self.quantities['xmomentum']
661            vh = self.quantities['ymomentum']
662
663            Vuh = uh.get_values(location='vertices', indices=[k])[0]
664            Euh = uh.get_values(location='edges', indices=[k])[0]
665            Cuh = uh.get_values(location='centroids', indices=[k])
666
667            Vvh = vh.get_values(location='vertices', indices=[k])[0]
668            Evh = vh.get_values(location='edges', indices=[k])[0]
669            Cvh = vh.get_values(location='centroids', indices=[k])
670
671            # Speeds in each direction
672            Vu = Vuh/(Vh + epsilon)
673            Eu = Euh/(Eh + epsilon)
674            Cu = Cuh/(Ch + epsilon)
675            name = 'U'
676            message  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n' \
677                 % (name.ljust(qwidth), Vu[0], Vu[1], Vu[2])
678
679            message += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n' \
680                 % (name.ljust(qwidth), Eu[0], Eu[1], Eu[2])
681
682            message += '    %s: centroid_value = %.4f\n' \
683                 % (name.ljust(qwidth), Cu[0])
684
685            msg += message
686
687            Vv = Vvh/(Vh + epsilon)
688            Ev = Evh/(Eh + epsilon)
689            Cv = Cvh/(Ch + epsilon)
690            name = 'V'
691            message  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n' \
692                 % (name.ljust(qwidth), Vv[0], Vv[1], Vv[2])
693
694            message += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n' \
695                 % (name.ljust(qwidth), Ev[0], Ev[1], Ev[2])
696
697            message += '    %s: centroid_value = %.4f\n'\
698                 %(name.ljust(qwidth), Cv[0])
699
700            msg += message
701
702            # Froude number in each direction
703            name = 'Froude (x)'
704            Vfx = Vu/(num.sqrt(g*Vh) + epsilon)
705            Efx = Eu/(num.sqrt(g*Eh) + epsilon)
706            Cfx = Cu/(num.sqrt(g*Ch) + epsilon)
707
708            message  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
709                 % (name.ljust(qwidth), Vfx[0], Vfx[1], Vfx[2])
710
711            message += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
712                 % (name.ljust(qwidth), Efx[0], Efx[1], Efx[2])
713
714            message += '    %s: centroid_value = %.4f\n'\
715                 % (name.ljust(qwidth), Cfx[0])
716
717            msg += message
718
719            name = 'Froude (y)'
720            Vfy = Vv/(num.sqrt(g*Vh) + epsilon)
721            Efy = Ev/(num.sqrt(g*Eh) + epsilon)
722            Cfy = Cv/(num.sqrt(g*Ch) + epsilon)
723
724            message  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
725                 % (name.ljust(qwidth), Vfy[0], Vfy[1], Vfy[2])
726
727            message += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
728                 % (name.ljust(qwidth), Efy[0], Efy[1], Efy[2])
729
730            message += '    %s: centroid_value = %.4f\n'\
731                 % (name.ljust(qwidth), Cfy[0])
732
733            msg += message
734
735        return msg
736       
737       
738
739    def compute_boundary_flows(self):
740        """Compute boundary flows at current timestep.
741                       
742        Quantities computed are:
743           Total inflow across boundary
744           Total outflow across boundary
745           Flow across each tagged boundary segment
746        """
747               
748        # Run through boundary array and compute for each segment
749        # the normal momentum ((uh, vh) dot normal) times segment length.
750        # Based on sign accumulate this into boundary_inflow and
751        # boundary_outflow.
752                       
753        # Compute flows along boundary
754       
755        uh = self.get_quantity('xmomentum').get_values(location='edges')
756        vh = self.get_quantity('ymomentum').get_values(location='edges')       
757       
758        # Loop through edges that lie on the boundary and calculate
759        # flows
760        boundary_flows = {}
761        total_boundary_inflow = 0.0
762        total_boundary_outflow = 0.0
763        for vol_id, edge_id in self.boundary:
764            # Compute normal flow across edge. Since normal vector points
765            # away from triangle, a positive sign means that water
766            # flows *out* from this triangle.
767             
768            momentum = [uh[vol_id, edge_id], vh[vol_id, edge_id]]
769            normal = self.mesh.get_normal(vol_id, edge_id)
770            length = self.mesh.get_edgelength(vol_id, edge_id)           
771            normal_flow = num.dot(momentum, normal)*length
772           
773            # Reverse sign so that + is taken to mean inflow
774            # and - means outflow. This is more intuitive.
775            edge_flow = -normal_flow
776           
777            # Tally up inflows and outflows separately
778            if edge_flow > 0:
779                # Flow is inflow     
780                total_boundary_inflow += edge_flow       
781            else:
782                # Flow is outflow
783                total_boundary_outflow += edge_flow   
784
785            # Tally up flows by boundary tag
786            tag = self.boundary[(vol_id, edge_id)]
787           
788            if tag not in boundary_flows:
789                boundary_flows[tag] = 0.0
790            boundary_flows[tag] += edge_flow
791           
792               
793        return boundary_flows, total_boundary_inflow, total_boundary_outflow
794       
795
796    def compute_forcing_flows(self):
797        """Compute flows in and out of domain due to forcing terms.
798                       
799        Quantities computed are:
800               
801       
802           Total inflow through forcing terms
803           Total outflow through forcing terms
804           Current total volume in domain       
805
806        """
807
808        #FIXME(Ole): We need to separate what part of explicit_update was
809        # due to the normal flux calculations and what is due to forcing terms.
810       
811        pass
812                       
813       
814    def compute_total_volume(self):
815        """Compute total volume (m^3) of water in entire domain
816        """
817       
818        area = self.mesh.get_areas()
819       
820        stage = self.get_quantity('stage').get_values(location='centroids')
821        elevation = \
822            self.get_quantity('elevation').get_values(location='centroids')
823        depth = stage-elevation
824       
825        return num.sum(depth*area)
826       
827       
828    def volumetric_balance_statistics(self):               
829        """Create volumetric balance report suitable for printing or logging.
830        """
831       
832        (boundary_flows, total_boundary_inflow,
833         total_boundary_outflow) = self.compute_boundary_flows() 
834       
835        message = '---------------------------\n'       
836        message += 'Volumetric balance report:\n'
837        message += '--------------------------\n'
838        message += 'Total boundary inflow [m^3/s]: %.2f\n' % total_boundary_inflow
839        message += 'Total boundary outflow [m^3/s]: %.2f\n' % total_boundary_outflow
840        message += 'Net boundary flow by tags [m^3/s]\n'
841        for tag in boundary_flows:
842            message += '    %s [m^3/s]: %.2f\n' % (tag, boundary_flows[tag])
843       
844        message += 'Total net boundary flow [m^3/s]: %.2f\n' % \
845                    (total_boundary_inflow + total_boundary_outflow) 
846        message += 'Total volume in domain [m^3]: %.2f\n' % \
847                    self.compute_total_volume()
848       
849        # The go through explicit forcing update and record the rate of change
850        # for stage and
851        # record into forcing_inflow and forcing_outflow. Finally compute 
852        # integral of depth to obtain total volume of domain.
853       
854        # FIXME(Ole): This part is not yet done.               
855       
856        return message       
857           
858################################################################################
859# End of class Shallow Water Domain
860################################################################################
861
862#-----------------
863# Flux computation
864#-----------------
865
866## @brief Compute fluxes and timestep suitable for all volumes in domain.
867# @param domain The domain to calculate fluxes for.
868def compute_fluxes(domain):
869    """Compute fluxes and timestep suitable for all volumes in domain.
870
871    Compute total flux for each conserved quantity using "flux_function"
872
873    Fluxes across each edge are scaled by edgelengths and summed up
874    Resulting flux is then scaled by area and stored in
875    explicit_update for each of the three conserved quantities
876    stage, xmomentum and ymomentum
877
878    The maximal allowable speed computed by the flux_function for each volume
879    is converted to a timestep that must not be exceeded. The minimum of
880    those is computed as the next overall timestep.
881
882    Post conditions:
883      domain.explicit_update is reset to computed flux values
884      domain.timestep is set to the largest step satisfying all volumes.
885
886    This wrapper calls the underlying C version of compute fluxes
887    """
888
889    import sys
890    from shallow_water_ext import compute_fluxes_ext_central \
891                                  as compute_fluxes_ext
892
893    # Shortcuts
894    Stage = domain.quantities['stage']
895    Xmom = domain.quantities['xmomentum']
896    Ymom = domain.quantities['ymomentum']
897    Bed = domain.quantities['elevation']
898
899    timestep = float(sys.maxint)
900
901    flux_timestep = compute_fluxes_ext(timestep,
902                                       domain.epsilon,
903                                       domain.H0,
904                                       domain.g,
905                                       domain.neighbours,
906                                       domain.neighbour_edges,
907                                       domain.normals,
908                                       domain.edgelengths,
909                                       domain.radii,
910                                       domain.areas,
911                                       domain.tri_full_flag,
912                                       Stage.edge_values,
913                                       Xmom.edge_values,
914                                       Ymom.edge_values,
915                                       Bed.edge_values,
916                                       Stage.boundary_values,
917                                       Xmom.boundary_values,
918                                       Ymom.boundary_values,
919                                       Stage.explicit_update,
920                                       Xmom.explicit_update,
921                                       Ymom.explicit_update,
922                                       domain.already_computed_flux,
923                                       domain.max_speed,
924                                       int(domain.optimise_dry_cells))
925
926    domain.flux_timestep = flux_timestep
927
928################################################################################
929# Module functions for gradient limiting
930################################################################################
931
932##
933# @brief Wrapper for C version of extrapolate_second_order_sw.
934# @param domain The domain to operate on.
935# @note MH090605 The following method belongs to the shallow_water domain class
936#       see comments in the corresponding method in shallow_water_ext.c
937def extrapolate_second_order_sw(domain):
938    """Wrapper calling C version of extrapolate_second_order_sw"""
939
940    from shallow_water_ext import extrapolate_second_order_sw as extrapol2
941
942    # Shortcuts
943    Stage = domain.quantities['stage']
944    Xmom = domain.quantities['xmomentum']
945    Ymom = domain.quantities['ymomentum']
946    Elevation = domain.quantities['elevation']
947
948    extrapol2(domain,
949              domain.surrogate_neighbours,
950              domain.number_of_boundaries,
951              domain.centroid_coordinates,
952              Stage.centroid_values,
953              Xmom.centroid_values,
954              Ymom.centroid_values,
955              Elevation.centroid_values,
956              domain.vertex_coordinates,
957              Stage.vertex_values,
958              Xmom.vertex_values,
959              Ymom.vertex_values,
960              Elevation.vertex_values,
961              int(domain.optimise_dry_cells))
962
963##
964# @brief Distribution from centroids to vertices specific to the SWW eqn.
965# @param domain The domain to operate on.
966def distribute_using_vertex_limiter(domain):
967    """Distribution from centroids to vertices specific to the SWW equation.
968
969    It will ensure that h (w-z) is always non-negative even in the
970    presence of steep bed-slopes by taking a weighted average between shallow
971    and deep cases.
972
973    In addition, all conserved quantities get distributed as per either a
974    constant (order==1) or a piecewise linear function (order==2).
975
976    FIXME: more explanation about removal of artificial variability etc
977
978    Precondition:
979      All quantities defined at centroids and bed elevation defined at
980      vertices.
981
982    Postcondition
983      Conserved quantities defined at vertices
984    """
985
986    # Remove very thin layers of water
987    protect_against_infinitesimal_and_negative_heights(domain)
988
989    # Extrapolate all conserved quantities
990    if domain.optimised_gradient_limiter:
991        # MH090605 if second order,
992        # perform the extrapolation and limiting on
993        # all of the conserved quantities
994
995        if (domain._order_ == 1):
996            for name in domain.conserved_quantities:
997                Q = domain.quantities[name]
998                Q.extrapolate_first_order()
999        elif domain._order_ == 2:
1000            domain.extrapolate_second_order_sw()
1001        else:
1002            raise Exception('Unknown order')
1003    else:
1004        # Old code:
1005        for name in domain.conserved_quantities:
1006            Q = domain.quantities[name]
1007
1008            if domain._order_ == 1:
1009                Q.extrapolate_first_order()
1010            elif domain._order_ == 2:
1011                Q.extrapolate_second_order_and_limit_by_vertex()
1012            else:
1013                raise Exception('Unknown order')
1014
1015    # Take bed elevation into account when water heights are small
1016    balance_deep_and_shallow(domain)
1017
1018    # Compute edge values by interpolation
1019    for name in domain.conserved_quantities:
1020        Q = domain.quantities[name]
1021        Q.interpolate_from_vertices_to_edges()
1022
1023##
1024# @brief Distribution from centroids to edges specific to the SWW eqn.
1025# @param domain The domain to operate on.
1026def distribute_using_edge_limiter(domain):
1027    """Distribution from centroids to edges specific to the SWW eqn.
1028
1029    It will ensure that h (w-z) is always non-negative even in the
1030    presence of steep bed-slopes by taking a weighted average between shallow
1031    and deep cases.
1032
1033    In addition, all conserved quantities get distributed as per either a
1034    constant (order==1) or a piecewise linear function (order==2).
1035
1036
1037    Precondition:
1038      All quantities defined at centroids and bed elevation defined at
1039      vertices.
1040
1041    Postcondition
1042      Conserved quantities defined at vertices
1043    """
1044
1045    # Remove very thin layers of water
1046    protect_against_infinitesimal_and_negative_heights(domain)
1047
1048    for name in domain.conserved_quantities:
1049        Q = domain.quantities[name]
1050        if domain._order_ == 1:
1051            Q.extrapolate_first_order()
1052        elif domain._order_ == 2:
1053            Q.extrapolate_second_order_and_limit_by_edge()
1054        else:
1055            raise Exception('Unknown order')
1056
1057    balance_deep_and_shallow(domain)
1058
1059    # Compute edge values by interpolation
1060    for name in domain.conserved_quantities:
1061        Q = domain.quantities[name]
1062        Q.interpolate_from_vertices_to_edges()
1063
1064##
1065# @brief  Protect against infinitesimal heights and associated high velocities.
1066# @param domain The domain to operate on.
1067def protect_against_infinitesimal_and_negative_heights(domain):
1068    """Protect against infinitesimal heights and associated high velocities"""
1069
1070    from shallow_water_ext import protect
1071
1072    # Shortcuts
1073    wc = domain.quantities['stage'].centroid_values
1074    zc = domain.quantities['elevation'].centroid_values
1075    xmomc = domain.quantities['xmomentum'].centroid_values
1076    ymomc = domain.quantities['ymomentum'].centroid_values
1077
1078    protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
1079            domain.epsilon, wc, zc, xmomc, ymomc)
1080
1081##
1082# @brief Wrapper for C function balance_deep_and_shallow_c().
1083# @param domain The domain to operate on.
1084def balance_deep_and_shallow(domain):
1085    """Compute linear combination between stage as computed by
1086    gradient-limiters limiting using w, and stage computed by
1087    gradient-limiters limiting using h (h-limiter).
1088    The former takes precedence when heights are large compared to the
1089    bed slope while the latter takes precedence when heights are
1090    relatively small.  Anything in between is computed as a balanced
1091    linear combination in order to avoid numerical disturbances which
1092    would otherwise appear as a result of hard switching between
1093    modes.
1094
1095    Wrapper for C implementation
1096    """
1097
1098    from shallow_water_ext import balance_deep_and_shallow \
1099                                  as balance_deep_and_shallow_c
1100
1101    # Shortcuts
1102    wc = domain.quantities['stage'].centroid_values
1103    zc = domain.quantities['elevation'].centroid_values
1104    wv = domain.quantities['stage'].vertex_values
1105    zv = domain.quantities['elevation'].vertex_values
1106
1107    # Momentums at centroids
1108    xmomc = domain.quantities['xmomentum'].centroid_values
1109    ymomc = domain.quantities['ymomentum'].centroid_values
1110
1111    # Momentums at vertices
1112    xmomv = domain.quantities['xmomentum'].vertex_values
1113    ymomv = domain.quantities['ymomentum'].vertex_values
1114
1115    balance_deep_and_shallow_c(domain,
1116                               wc, zc, wv, zv, wc,
1117                               xmomc, ymomc, xmomv, ymomv)
1118
1119
1120
1121################################################################################
1122# Standard forcing terms
1123################################################################################
1124
1125##
1126# @brief Apply gravitational pull in the presence of bed slope.
1127# @param domain The domain to apply gravity to.
1128# @note Wrapper for C function gravity_c().
1129def gravity(domain):
1130    """Apply gravitational pull in the presence of bed slope
1131    Wrapper calls underlying C implementation
1132    """
1133
1134    from shallow_water_ext import gravity as gravity_c
1135
1136    xmom_update = domain.quantities['xmomentum'].explicit_update
1137    ymom_update = domain.quantities['ymomentum'].explicit_update
1138
1139    stage = domain.quantities['stage']
1140    elevation = domain.quantities['elevation']
1141
1142    #FIXME SR Should avoid allocating memory!
1143    height = stage.centroid_values - elevation.centroid_values
1144    elevation = elevation.vertex_values
1145
1146    point = domain.get_vertex_coordinates()
1147
1148    gravity_c(domain.g, height, elevation, point, xmom_update, ymom_update)
1149
1150##
1151# @brief Apply friction to a surface (implicit).
1152# @param domain The domain to apply Manning friction to.
1153# @note Wrapper for C function manning_friction_c().
1154def manning_friction_implicit(domain):
1155    """Apply (Manning) friction to water momentum
1156    Wrapper for c version
1157    """
1158
1159    from shallow_water_ext import manning_friction_old
1160    from shallow_water_ext import manning_friction_new
1161
1162    xmom = domain.quantities['xmomentum']
1163    ymom = domain.quantities['ymomentum']
1164
1165    x = domain.get_vertex_coordinates()
1166   
1167    w = domain.quantities['stage'].centroid_values
1168    z = domain.quantities['elevation'].vertex_values
1169
1170    uh = xmom.centroid_values
1171    vh = ymom.centroid_values
1172    eta = domain.quantities['friction'].centroid_values
1173
1174    xmom_update = xmom.semi_implicit_update
1175    ymom_update = ymom.semi_implicit_update
1176
1177    eps = domain.minimum_allowed_height
1178    g = domain.g
1179
1180    if domain.use_new_mannings:
1181        manning_friction_new(g, eps, x, w, uh, vh, z, eta, xmom_update, \
1182                                ymom_update)
1183    else:
1184        manning_friction_old(g, eps, w, uh, vh, z, eta, xmom_update, \
1185                                ymom_update)
1186   
1187
1188##
1189# @brief Apply friction to a surface (explicit).
1190# @param domain The domain to apply Manning friction to.
1191# @note Wrapper for C function manning_friction_c().
1192def manning_friction_explicit(domain):
1193    """Apply (Manning) friction to water momentum
1194    Wrapper for c version
1195    """
1196
1197    from shallow_water_ext import manning_friction_old
1198    from shallow_water_ext import manning_friction_new
1199
1200    xmom = domain.quantities['xmomentum']
1201    ymom = domain.quantities['ymomentum']
1202
1203    x = domain.get_vertex_coordinates()
1204   
1205    w = domain.quantities['stage'].centroid_values
1206    z = domain.quantities['elevation'].vertex_values
1207
1208    uh = xmom.centroid_values
1209    vh = ymom.centroid_values
1210    eta = domain.quantities['friction'].centroid_values
1211
1212    xmom_update = xmom.explicit_update
1213    ymom_update = ymom.explicit_update
1214
1215    eps = domain.minimum_allowed_height
1216
1217    if domain.use_new_mannings:
1218        manning_friction_new(domain.g, eps, x, w, uh, vh, z, eta, xmom_update, \
1219                            ymom_update)
1220    else:
1221        manning_friction_old(domain.g, eps, w, uh, vh, z, eta, xmom_update, \
1222                            ymom_update)
1223
1224
1225
1226# FIXME (Ole): This was implemented for use with one of the analytical solutions
1227##
1228# @brief Apply linear friction to a surface.
1229# @param domain The domain to apply Manning friction to.
1230# @note Is this still used (30 Oct 2007)?
1231def linear_friction(domain):
1232    """Apply linear friction to water momentum
1233
1234    Assumes quantity: 'linear_friction' to be present
1235    """
1236
1237    w = domain.quantities['stage'].centroid_values
1238    z = domain.quantities['elevation'].centroid_values
1239    h = w-z
1240
1241    uh = domain.quantities['xmomentum'].centroid_values
1242    vh = domain.quantities['ymomentum'].centroid_values
1243    tau = domain.quantities['linear_friction'].centroid_values
1244
1245    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1246    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1247
1248    num_tris = len(domain)
1249    eps = domain.minimum_allowed_height
1250
1251    for k in range(num_tris):
1252        if tau[k] >= eps:
1253            if h[k] >= eps:
1254                S = -tau[k]/h[k]
1255
1256                #Update momentum
1257                xmom_update[k] += S*uh[k]
1258                ymom_update[k] += S*vh[k]
1259
1260def depth_dependent_friction(domain, default_friction,
1261                             surface_roughness_data,
1262                             verbose=False):
1263    """Returns an array of friction values for each wet element adjusted for
1264            depth.
1265
1266    Inputs:
1267        domain - computational domain object
1268        default_friction - depth independent bottom friction
1269        surface_roughness_data - N x 5 array of n0, d1, n1, d2, n2 values
1270        for each friction region.
1271
1272    Outputs:
1273        wet_friction - Array that can be used directly to update friction as
1274                        follows:
1275                       domain.set_quantity('friction', wet_friction)
1276
1277       
1278       
1279    """
1280   
1281    default_n0 = 0  # James - this was missing, don't know what it should be
1282   
1283    # Create a temp array to store updated depth dependent
1284    # friction for wet elements
1285    # EHR this is outwardly inneficient but not obvious how to avoid
1286    # recreating each call??????
1287
1288    wet_friction    = num.zeros(len(domain), num.float)
1289    wet_friction[:] = default_n0  # Initially assign default_n0 to all array so
1290                                  # sure have no zeros values
1291   
1292    # create depth instance for this timestep   
1293    depth = domain.create_quantity_from_expression('stage - elevation') 
1294    # Recompute depth as vector 
1295    d_vals = depth.get_values(location='centroids')
1296 
1297    # rebuild the 'friction' values adjusted for depth at this instant
1298    # loop for each wet element in domain
1299   
1300    for i in domain.get_wet_elements():       
1301        # Get roughness data for each element
1302        d1 = float(surface_roughness_data[i, 1])
1303        n1 = float(surface_roughness_data[i, 2])
1304        d2 = float(surface_roughness_data[i, 3])
1305        n2 = float(surface_roughness_data[i, 4])
1306       
1307       
1308        # Recompute friction values from depth for this element
1309               
1310        if d_vals[i]   <= d1: 
1311            ddf = n1
1312        elif d_vals[i] >= d2:
1313            ddf = n2
1314        else:
1315            ddf = n1 + ((n2-n1)/(d2-d1))*(d_vals[i]-d1)
1316           
1317        # check sanity of result
1318        if (ddf  < 0.010 or \
1319                            ddf > 9999.0) :
1320            log.critical('>>>> WARNING: computed depth_dependent friction '
1321                         'out of range, ddf%f, n1=%f, n2=%f'
1322                         % (ddf, n1, n2))
1323       
1324        # update depth dependent friction  for that wet element
1325        wet_friction[i] = ddf
1326       
1327    # EHR add code to show range of 'friction across domain at this instant as
1328    # sanity check?????????
1329   
1330    if verbose :
1331        # return array of domain nvals
1332        nvals = domain.get_quantity('friction').get_values(location='centroids')
1333        n_min = min(nvals)
1334        n_max = max(nvals)
1335       
1336        log.critical('         ++++ calculate_depth_dependent_friction - '
1337                     'Updated friction - range  %7.3f to %7.3f'
1338                     % (n_min, n_max))
1339   
1340    return wet_friction
1341
1342
1343
1344################################################################################
1345# Initialise module
1346################################################################################
1347
1348def _raise_compile_exception():
1349    """ Raise exception if compiler not available. """
1350    msg = 'C implementations could not be accessed by %s.\n ' % __file__
1351    msg += 'Make sure compile_all.py has been run as described in '
1352    msg += 'the ANUGA installation guide.'
1353    raise Exception(msg)
1354
1355from anuga.utilities import compile
1356if not compile.can_use_C_extension('shallow_water_ext.c'):
1357    _raise_compile_exception()
1358
1359if __name__ == "__main__":
1360    pass
Note: See TracBrowser for help on using the repository browser.