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

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

Changes to address ticket 360.

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