source: anuga_core/source/anuga/shallow_water/shallow_water_domain.py @ 5294

Last change on this file since 5294 was 5294, checked in by ole, 16 years ago

Work on Inflow and Rain forcing terms.
Also looked at a few backwards compatibility issues with a view to
commit new limiter options

  • Property svn:keywords set to LastChangedDate LastChangedRevision LastChangedBy HeadURL Id
File size: 58.8 KB
Line 
1"""Finite-volume computations of the shallow water wave equation.
2
3Title: ANGUA shallow_water_domain - 2D triangular domains for finite-volume
4       computations of the shallow water wave equation.
5
6
7Author: Ole Nielsen, Ole.Nielsen@ga.gov.au
8        Stephen Roberts, Stephen.Roberts@anu.edu.au
9        Duncan Gray, Duncan.Gray@ga.gov.au
10
11CreationDate: 2004
12
13Description:
14    This module contains a specialisation of class Domain from
15    module domain.py consisting of methods specific to the
16    Shallow Water Wave Equation
17
18    U_t + E_x + G_y = S
19
20    where
21
22    U = [w, uh, vh]
23    E = [uh, u^2h + gh^2/2, uvh]
24    G = [vh, uvh, v^2h + gh^2/2]
25    S represents source terms forcing the system
26    (e.g. gravity, friction, wind stress, ...)
27
28    and _t, _x, _y denote the derivative with respect to t, x and y
29    respectively.
30
31
32    The quantities are
33
34    symbol    variable name    explanation
35    x         x                horizontal distance from origin [m]
36    y         y                vertical distance from origin [m]
37    z         elevation        elevation of bed on which flow is modelled [m]
38    h         height           water height above z [m]
39    w         stage            absolute water level, w = z+h [m]
40    u                          speed in the x direction [m/s]
41    v                          speed in the y direction [m/s]
42    uh        xmomentum        momentum in the x direction [m^2/s]
43    vh        ymomentum        momentum in the y direction [m^2/s]
44
45    eta                        mannings friction coefficient [to appear]
46    nu                         wind stress coefficient [to appear]
47
48    The conserved quantities are w, uh, vh
49
50Reference:
51    Catastrophic Collapse of Water Supply Reservoirs in Urban Areas,
52    Christopher Zoppou and Stephen Roberts,
53    Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999
54
55    Hydrodynamic modelling of coastal inundation.
56    Nielsen, O., S. Roberts, D. Gray, A. McPherson and A. Hitchman
57    In Zerger, A. and Argent, R.M. (eds) MODSIM 2005 International Congress on
58    Modelling and Simulation. Modelling and Simulation Society of Australia and
59    New Zealand, December 2005, pp. 518-523. ISBN: 0-9758400-2-9.
60    http://www.mssanz.org.au/modsim05/papers/nielsen.pdf
61
62
63SeeAlso:
64    TRAC administration of ANUGA (User Manuals etc) at
65    https://datamining.anu.edu.au/anuga and Subversion repository at
66    $HeadURL: anuga_core/source/anuga/shallow_water/shallow_water_domain.py $
67
68Constraints: See GPL license in the user guide
69Version: 1.0 ($Revision: 5294 $)
70ModifiedBy:
71    $Author: ole $
72    $Date: 2008-05-08 12:58:15 +0000 (Thu, 08 May 2008) $
73
74"""
75
76# Subversion keywords:
77#
78# $LastChangedDate: 2008-05-08 12:58:15 +0000 (Thu, 08 May 2008) $
79# $LastChangedRevision: 5294 $
80# $LastChangedBy: ole $
81
82from Numeric import zeros, ones, Float, array, sum, size
83from Numeric import compress, arange
84
85
86from anuga.abstract_2d_finite_volumes.domain import Domain as Generic_Domain
87from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
88     import Boundary
89from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
90     import File_boundary
91from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
92     import Dirichlet_boundary
93from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
94     import Time_boundary
95from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
96     import Transmissive_boundary
97
98from anuga.utilities.numerical_tools import gradient, mean, ensure_numeric
99from anuga.config import minimum_storable_height
100from anuga.config import minimum_allowed_height, maximum_allowed_speed
101from anuga.config import g, epsilon, beta_h, beta_w, beta_w_dry,\
102     beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, tight_slope_limiters
103from anuga.config import alpha_balance
104from anuga.config import optimise_dry_cells
105from anuga.config import optimised_gradient_limiter
106from anuga.config import use_edge_limiter
107from anuga.config import use_centroid_velocities
108
109
110from anuga.utilities.polygon import inside_polygon, polygon_area       
111
112
113
114#---------------------
115# Shallow water domain
116#---------------------
117class Domain(Generic_Domain):
118
119    conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
120    other_quantities = ['elevation', 'friction']
121   
122    def __init__(self,
123                 coordinates=None,
124                 vertices=None,
125                 boundary=None,
126                 tagged_elements=None,
127                 geo_reference=None,
128                 use_inscribed_circle=False,
129                 mesh_filename=None,
130                 use_cache=False,
131                 verbose=False,
132                 full_send_dict=None,
133                 ghost_recv_dict=None,
134                 processor=0,
135                 numproc=1,
136                 number_of_full_nodes=None,
137                 number_of_full_triangles=None):
138
139
140        other_quantities = ['elevation', 'friction']
141        Generic_Domain.__init__(self,
142                                coordinates,
143                                vertices,
144                                boundary,
145                                Domain.conserved_quantities,
146                                Domain.other_quantities,
147                                tagged_elements,
148                                geo_reference,
149                                use_inscribed_circle,
150                                mesh_filename,
151                                use_cache,
152                                verbose,
153                                full_send_dict,
154                                ghost_recv_dict,
155                                processor,
156                                numproc,
157                                number_of_full_nodes=number_of_full_nodes,
158                                number_of_full_triangles=number_of_full_triangles) 
159
160
161        self.set_minimum_allowed_height(minimum_allowed_height)
162       
163        self.maximum_allowed_speed = maximum_allowed_speed
164        self.g = g
165        self.beta_w      = beta_w
166        self.beta_w_dry  = beta_w_dry
167        self.beta_uh     = beta_uh
168        self.beta_uh_dry = beta_uh_dry
169        self.beta_vh     = beta_vh
170        self.beta_vh_dry = beta_vh_dry
171        self.beta_h      = beta_h
172        self.alpha_balance = alpha_balance
173
174        self.tight_slope_limiters = tight_slope_limiters
175        self.optimise_dry_cells = optimise_dry_cells
176
177        self.forcing_terms.append(manning_friction_implicit)
178        self.forcing_terms.append(gravity)
179
180        # Stored output
181        self.store = True
182        self.format = 'sww'
183        self.set_store_vertices_uniquely(False)
184        self.minimum_storable_height = minimum_storable_height
185        self.quantities_to_be_stored = ['stage','xmomentum','ymomentum']
186
187        # Limiters
188        self.use_edge_limiter = use_edge_limiter
189        self.optimised_gradient_limiter = optimised_gradient_limiter
190        self.use_centroid_velocities = use_centroid_velocities
191
192
193    def set_all_limiters(self, beta):
194        """Shorthand to assign one constant value [0,1[ to all limiters.
195        0 Corresponds to first order, where as larger values make use of
196        the second order scheme.
197        """
198
199        self.beta_w      = beta
200        self.beta_w_dry  = beta
201        self.quantities['stage'].beta = beta
202       
203        self.beta_uh     = beta
204        self.beta_uh_dry = beta
205        self.quantities['xmomentum'].beta = beta
206       
207        self.beta_vh     = beta
208        self.beta_vh_dry = beta
209        self.quantities['ymomentum'].beta = beta
210       
211        self.beta_h      = beta
212       
213
214    def set_store_vertices_uniquely(self, flag, reduction=None):
215        """Decide whether vertex values should be stored uniquely as
216        computed in the model or whether they should be reduced to one
217        value per vertex using self.reduction.
218        """
219
220        # FIXME (Ole): how about using the word continuous vertex values?
221        self.smooth = not flag
222
223        # Reduction operation for get_vertex_values
224        if reduction is None:
225            self.reduction = mean
226            #self.reduction = min  #Looks better near steep slopes
227
228
229    def set_minimum_storable_height(self, minimum_storable_height):
230        """
231        Set the minimum depth that will be recognised when writing
232        to an sww file. This is useful for removing thin water layers
233        that seems to be caused by friction creep.
234
235        The minimum allowed sww depth is in meters.
236        """
237        self.minimum_storable_height = minimum_storable_height
238
239
240    def set_minimum_allowed_height(self, minimum_allowed_height):
241        """
242        Set the minimum depth that will be recognised in the numerical
243        scheme
244
245        The minimum allowed depth is in meters.
246
247        The parameter H0 (Minimal height for flux computation)
248        is also set by this function
249        """
250
251        #FIXME (Ole): rename H0 to minimum_allowed_height_in_flux_computation
252
253        #FIXME (Ole): Maybe use histogram to identify isolated extreme speeds
254        #and deal with them adaptively similarly to how we used to use 1 order
255        #steps to recover.
256        self.minimum_allowed_height = minimum_allowed_height
257        self.H0 = minimum_allowed_height   
258       
259
260    def set_maximum_allowed_speed(self, maximum_allowed_speed):
261        """
262        Set the maximum particle speed that is allowed in water
263        shallower than minimum_allowed_height. This is useful for
264        controlling speeds in very thin layers of water and at the same time
265        allow some movement avoiding pooling of water.
266
267        """
268        self.maximum_allowed_speed = maximum_allowed_speed
269
270
271    def set_points_file_block_line_size(self,points_file_block_line_size):
272        """
273        Set the minimum depth that will be recognised when writing
274        to an sww file. This is useful for removing thin water layers
275        that seems to be caused by friction creep.
276
277        The minimum allowed sww depth is in meters.
278        """
279        self.points_file_block_line_size = points_file_block_line_size
280       
281       
282    def set_quantities_to_be_stored(self, q):
283        """Specify which quantities will be stored in the sww file.
284
285        q must be either:
286          - the name of a quantity
287          - a list of quantity names
288          - None
289
290        In the two first cases, the named quantities will be stored at
291        each yieldstep (This is in addition to the quantities elevation
292        and friction)
293       
294        If q is None, storage will be switched off altogether.
295        """
296
297
298        if q is None:
299            self.quantities_to_be_stored = []
300            self.store = False
301            return
302
303        if isinstance(q, basestring):
304            q = [q] # Turn argument into a list
305
306        # Check correcness
307        for quantity_name in q:
308            msg = 'Quantity %s is not a valid conserved quantity'\
309                  %quantity_name
310           
311            assert quantity_name in self.conserved_quantities, msg
312
313        self.quantities_to_be_stored = q
314
315
316
317    def get_wet_elements(self, indices=None):
318        """Return indices for elements where h > minimum_allowed_height
319
320        Optional argument:
321            indices is the set of element ids that the operation applies to.
322
323        Usage:
324            indices = get_wet_elements()
325
326        Note, centroid values are used for this operation           
327        """
328
329        # Water depth below which it is considered to be 0 in the model
330        # FIXME (Ole): Allow this to be specified as a keyword argument as well
331        from anuga.config import minimum_allowed_height
332
333       
334        elevation = self.get_quantity('elevation').\
335                    get_values(location='centroids', indices=indices)
336        stage = self.get_quantity('stage').\
337                get_values(location='centroids', indices=indices)       
338        depth = stage - elevation
339
340        # Select indices for which depth > 0
341        wet_indices = compress(depth > minimum_allowed_height,
342                               arange(len(depth)))
343        return wet_indices
344
345
346    def get_maximum_inundation_elevation(self, indices=None):
347        """Return highest elevation where h > 0
348
349        Optional argument:
350            indices is the set of element ids that the operation applies to.
351
352        Usage:
353            q = get_maximum_inundation_elevation()
354
355        Note, centroid values are used for this operation           
356        """
357
358        wet_elements = self.get_wet_elements(indices)
359        return self.get_quantity('elevation').\
360               get_maximum_value(indices=wet_elements)
361
362
363    def get_maximum_inundation_location(self, indices=None):
364        """Return location of highest elevation where h > 0
365
366        Optional argument:
367            indices is the set of element ids that the operation applies to.
368
369        Usage:
370            q = get_maximum_inundation_location()
371
372        Note, centroid values are used for this operation           
373        """
374
375        wet_elements = self.get_wet_elements(indices)
376        return self.get_quantity('elevation').\
377               get_maximum_location(indices=wet_elements)   
378
379    def check_integrity(self):
380        Generic_Domain.check_integrity(self)
381
382        #Check that we are solving the shallow water wave equation
383
384        msg = 'First conserved quantity must be "stage"'
385        assert self.conserved_quantities[0] == 'stage', msg
386        msg = 'Second conserved quantity must be "xmomentum"'
387        assert self.conserved_quantities[1] == 'xmomentum', msg
388        msg = 'Third conserved quantity must be "ymomentum"'
389        assert self.conserved_quantities[2] == 'ymomentum', msg
390
391    def extrapolate_second_order_sw(self):
392        #Call correct module function
393        #(either from this module or C-extension)
394        extrapolate_second_order_sw(self)
395
396    def compute_fluxes(self):
397        #Call correct module function
398        #(either from this module or C-extension)
399        compute_fluxes(self)
400
401    def distribute_to_vertices_and_edges(self):
402        # Call correct module function
403        # (either from this module or C-extension)
404        if self.use_edge_limiter:
405            protect_against_infinitesimal_and_negative_heights(self)
406            for name in self.conserved_quantities:
407                Q = self.quantities[name]
408                if self._order_ == 1:
409                    Q.extrapolate_first_order()
410                elif self._order_ == 2:
411                    Q.extrapolate_second_order_and_limit()
412                else:
413                    raise 'Unknown order'
414            balance_deep_and_shallow(self)
415
416            #Compute edge values by interpolation
417            for name in self.conserved_quantities:
418                Q = self.quantities[name]
419                Q.interpolate_from_vertices_to_edges()
420        else:
421            distribute_to_vertices_and_edges(self)
422
423
424
425
426    def evolve(self,
427               yieldstep = None,
428               finaltime = None,
429               duration = None,
430               skip_initial_step = False):
431        """Specialisation of basic evolve method from parent class
432        """
433
434        # Call check integrity here rather than from user scripts
435        # self.check_integrity()
436
437        msg = 'Parameter beta_h must be in the interval [0, 2['
438        assert 0 <= self.beta_h <= 2.0, msg
439        msg = 'Parameter beta_w must be in the interval [0, 2['
440        assert 0 <= self.beta_w <= 2.0, msg
441
442
443        # Initial update of vertex and edge values before any STORAGE
444        # and or visualisation
445        # This is done again in the initialisation of the Generic_Domain
446        # evolve loop but we do it here to ensure the values are ok for storage
447        self.distribute_to_vertices_and_edges()
448
449        if self.store is True and self.time == 0.0:
450            self.initialise_storage()
451            # print 'Storing results in ' + self.writer.filename
452        else:
453            pass
454            # print 'Results will not be stored.'
455            # print 'To store results set domain.store = True'
456            # FIXME: Diagnostic output should be controlled by
457            # a 'verbose' flag living in domain (or in a parent class)
458
459        # Call basic machinery from parent class
460        for t in Generic_Domain.evolve(self,
461                                       yieldstep=yieldstep,
462                                       finaltime=finaltime,
463                                       duration=duration,
464                                       skip_initial_step=skip_initial_step):
465
466            # Store model data, e.g. for subsequent visualisation
467            if self.store is True:
468                self.store_timestep(self.quantities_to_be_stored)
469
470            # FIXME: Could maybe be taken from specified list
471            # of 'store every step' quantities
472
473            # Pass control on to outer loop for more specific actions
474            yield(t)
475
476
477    def initialise_storage(self):
478        """Create and initialise self.writer object for storing data.
479        Also, save x,y and bed elevation
480        """
481
482        from anuga.shallow_water.data_manager import get_dataobject
483
484        # Initialise writer
485        self.writer = get_dataobject(self, mode = 'w')
486
487        # Store vertices and connectivity
488        self.writer.store_connectivity()
489
490
491    def store_timestep(self, name):
492        """Store named quantity and time.
493
494        Precondition:
495           self.write has been initialised
496        """
497        self.writer.store_timestep(name)
498
499       
500    def timestepping_statistics(self,
501                                track_speeds=False,
502                                triangle_id=None):       
503        """Return string with time stepping statistics for printing or logging
504
505        Optional boolean keyword track_speeds decides whether to report
506        location of smallest timestep as well as a histogram and percentile
507        report.
508        """
509
510        from Numeric import sqrt
511        from anuga.config import epsilon, g               
512
513
514        # Call basic machinery from parent class
515        msg = Generic_Domain.timestepping_statistics(self,
516                                                     track_speeds,
517                                                     triangle_id)
518
519        if track_speeds is True:
520
521            # qwidth determines the text field used for quantities
522            qwidth = self.qwidth
523       
524            # Selected triangle
525            k = self.k
526
527            # Report some derived quantities at vertices, edges and centroid
528            # specific to the shallow water wave equation
529
530            z = self.quantities['elevation']
531            w = self.quantities['stage']           
532
533            Vw = w.get_values(location='vertices', indices=[k])[0]
534            Ew = w.get_values(location='edges', indices=[k])[0]
535            Cw = w.get_values(location='centroids', indices=[k])
536
537            Vz = z.get_values(location='vertices', indices=[k])[0]
538            Ez = z.get_values(location='edges', indices=[k])[0]
539            Cz = z.get_values(location='centroids', indices=[k])                             
540               
541
542            name = 'depth'
543            Vh = Vw-Vz
544            Eh = Ew-Ez
545            Ch = Cw-Cz
546           
547            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
548                 %(name.ljust(qwidth), Vh[0], Vh[1], Vh[2])
549           
550            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
551                 %(name.ljust(qwidth), Eh[0], Eh[1], Eh[2])
552           
553            s += '    %s: centroid_value = %.4f\n'\
554                 %(name.ljust(qwidth), Ch[0])                               
555           
556            msg += s
557
558            uh = self.quantities['xmomentum']
559            vh = self.quantities['ymomentum']
560
561            Vuh = uh.get_values(location='vertices', indices=[k])[0]
562            Euh = uh.get_values(location='edges', indices=[k])[0]
563            Cuh = uh.get_values(location='centroids', indices=[k])
564           
565            Vvh = vh.get_values(location='vertices', indices=[k])[0]
566            Evh = vh.get_values(location='edges', indices=[k])[0]
567            Cvh = vh.get_values(location='centroids', indices=[k])
568
569            # Speeds in each direction
570            Vu = Vuh/(Vh + epsilon)
571            Eu = Euh/(Eh + epsilon)
572            Cu = Cuh/(Ch + epsilon)           
573            name = 'U'
574            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
575                 %(name.ljust(qwidth), Vu[0], Vu[1], Vu[2])
576           
577            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
578                 %(name.ljust(qwidth), Eu[0], Eu[1], Eu[2])
579           
580            s += '    %s: centroid_value = %.4f\n'\
581                 %(name.ljust(qwidth), Cu[0])                               
582           
583            msg += s
584
585            Vv = Vvh/(Vh + epsilon)
586            Ev = Evh/(Eh + epsilon)
587            Cv = Cvh/(Ch + epsilon)           
588            name = 'V'
589            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
590                 %(name.ljust(qwidth), Vv[0], Vv[1], Vv[2])
591           
592            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
593                 %(name.ljust(qwidth), Ev[0], Ev[1], Ev[2])
594           
595            s += '    %s: centroid_value = %.4f\n'\
596                 %(name.ljust(qwidth), Cv[0])                               
597           
598            msg += s
599
600
601            # Froude number in each direction
602            name = 'Froude (x)'
603            Vfx = Vu/(sqrt(g*Vh) + epsilon)
604            Efx = Eu/(sqrt(g*Eh) + epsilon)
605            Cfx = Cu/(sqrt(g*Ch) + epsilon)
606           
607            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
608                 %(name.ljust(qwidth), Vfx[0], Vfx[1], Vfx[2])
609           
610            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
611                 %(name.ljust(qwidth), Efx[0], Efx[1], Efx[2])
612           
613            s += '    %s: centroid_value = %.4f\n'\
614                 %(name.ljust(qwidth), Cfx[0])                               
615           
616            msg += s
617
618
619            name = 'Froude (y)'
620            Vfy = Vv/(sqrt(g*Vh) + epsilon)
621            Efy = Ev/(sqrt(g*Eh) + epsilon)
622            Cfy = Cv/(sqrt(g*Ch) + epsilon)
623           
624            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
625                 %(name.ljust(qwidth), Vfy[0], Vfy[1], Vfy[2])
626           
627            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
628                 %(name.ljust(qwidth), Efy[0], Efy[1], Efy[2])
629           
630            s += '    %s: centroid_value = %.4f\n'\
631                 %(name.ljust(qwidth), Cfy[0])                               
632           
633            msg += s           
634
635               
636
637        return msg
638       
639       
640
641#=============== End of class Shallow Water Domain ===============================
642
643
644#-----------------
645# Flux computation
646#-----------------
647
648def compute_fluxes(domain):
649    """Compute all fluxes and the timestep suitable for all volumes
650    in domain.
651
652    Compute total flux for each conserved quantity using "flux_function"
653
654    Fluxes across each edge are scaled by edgelengths and summed up
655    Resulting flux is then scaled by area and stored in
656    explicit_update for each of the three conserved quantities
657    stage, xmomentum and ymomentum
658
659    The maximal allowable speed computed by the flux_function for each volume
660    is converted to a timestep that must not be exceeded. The minimum of
661    those is computed as the next overall timestep.
662
663    Post conditions:
664      domain.explicit_update is reset to computed flux values
665      domain.timestep is set to the largest step satisfying all volumes.
666   
667
668    This wrapper calls the underlying C version of compute fluxes
669    """
670
671    import sys
672
673    N = len(domain) # number_of_triangles
674
675    # Shortcuts
676    Stage = domain.quantities['stage']
677    Xmom = domain.quantities['xmomentum']
678    Ymom = domain.quantities['ymomentum']
679    Bed = domain.quantities['elevation']
680
681    timestep = float(sys.maxint)
682    from shallow_water_ext import\
683         compute_fluxes_ext_central as compute_fluxes_ext
684
685
686    flux_timestep = compute_fluxes_ext(timestep,
687                                       domain.epsilon,
688                                       domain.H0,
689                                       domain.g,
690                                       domain.neighbours,
691                                       domain.neighbour_edges,
692                                       domain.normals,
693                                       domain.edgelengths,
694                                       domain.radii,
695                                       domain.areas,
696                                       domain.tri_full_flag,
697                                       Stage.edge_values,
698                                       Xmom.edge_values,
699                                       Ymom.edge_values,
700                                       Bed.edge_values,
701                                       Stage.boundary_values,
702                                       Xmom.boundary_values,
703                                       Ymom.boundary_values,
704                                       Stage.explicit_update,
705                                       Xmom.explicit_update,
706                                       Ymom.explicit_update,
707                                       domain.already_computed_flux,
708                                       domain.max_speed,
709                                       int(domain.optimise_dry_cells))
710
711    domain.flux_timestep = flux_timestep
712
713
714
715#---------------------------------------
716# Module functions for gradient limiting
717#---------------------------------------
718
719
720# MH090605 The following method belongs to the shallow_water domain class
721# see comments in the corresponding method in shallow_water_ext.c
722def extrapolate_second_order_sw(domain):
723    """Wrapper calling C version of extrapolate_second_order_sw
724    """
725    import sys
726
727    N = len(domain) # number_of_triangles
728
729    # Shortcuts
730    Stage = domain.quantities['stage']
731    Xmom = domain.quantities['xmomentum']
732    Ymom = domain.quantities['ymomentum']
733    Elevation = domain.quantities['elevation']
734
735    from shallow_water_ext import extrapolate_second_order_sw as extrapol2
736    extrapol2(domain,
737              domain.surrogate_neighbours,
738              domain.number_of_boundaries,
739              domain.centroid_coordinates,
740              Stage.centroid_values,
741              Xmom.centroid_values,
742              Ymom.centroid_values,
743              Elevation.centroid_values,
744              domain.vertex_coordinates,
745              Stage.vertex_values,
746              Xmom.vertex_values,
747              Ymom.vertex_values,
748              Elevation.vertex_values,
749              int(domain.optimise_dry_cells),
750              int(domain.use_centroid_velocities))
751
752
753def distribute_to_vertices_and_edges(domain):
754    """Distribution from centroids to vertices specific to the
755    shallow water wave
756    equation.
757
758    It will ensure that h (w-z) is always non-negative even in the
759    presence of steep bed-slopes by taking a weighted average between shallow
760    and deep cases.
761
762    In addition, all conserved quantities get distributed as per either a
763    constant (order==1) or a piecewise linear function (order==2).
764
765    FIXME: more explanation about removal of artificial variability etc
766
767    Precondition:
768      All quantities defined at centroids and bed elevation defined at
769      vertices.
770
771    Postcondition
772      Conserved quantities defined at vertices
773
774    """
775
776   
777
778    # Remove very thin layers of water
779    protect_against_infinitesimal_and_negative_heights(domain)
780
781    # Extrapolate all conserved quantities
782    if domain.optimised_gradient_limiter:
783        # MH090605 if second order,
784        # perform the extrapolation and limiting on
785        # all of the conserved quantities
786
787        if (domain._order_ == 1):
788            for name in domain.conserved_quantities:
789                Q = domain.quantities[name]
790                Q.extrapolate_first_order()
791        elif domain._order_ == 2:
792            domain.extrapolate_second_order_sw()
793        else:
794            raise 'Unknown order'
795    else:
796        # Old code:
797        for name in domain.conserved_quantities:
798            Q = domain.quantities[name]
799
800            if domain._order_ == 1:
801                Q.extrapolate_first_order()
802            elif domain._order_ == 2:
803                Q.extrapolate_second_order()
804                Q.limit()
805            else:
806                raise 'Unknown order'
807
808
809    # Take bed elevation into account when water heights are small
810    balance_deep_and_shallow(domain)
811
812    # Compute edge values by interpolation
813    for name in domain.conserved_quantities:
814        Q = domain.quantities[name]
815        Q.interpolate_from_vertices_to_edges()
816
817
818def protect_against_infinitesimal_and_negative_heights(domain):
819    """Protect against infinitesimal heights and associated high velocities
820    """
821
822    # Shortcuts
823    wc = domain.quantities['stage'].centroid_values
824    zc = domain.quantities['elevation'].centroid_values
825    xmomc = domain.quantities['xmomentum'].centroid_values
826    ymomc = domain.quantities['ymomentum'].centroid_values
827
828    from shallow_water_ext import protect
829
830    protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
831            domain.epsilon, wc, zc, xmomc, ymomc)
832
833
834def h_limiter(domain):
835    """Limit slopes for each volume to eliminate artificial variance
836    introduced by e.g. second order extrapolator
837
838    limit on h = w-z
839
840    This limiter depends on two quantities (w,z) so it resides within
841    this module rather than within quantity.py
842
843    Wrapper for c-extension
844    """
845
846    N = len(domain) # number_of_triangles
847    beta_h = domain.beta_h
848
849    # Shortcuts
850    wc = domain.quantities['stage'].centroid_values
851    zc = domain.quantities['elevation'].centroid_values
852    hc = wc - zc
853
854    wv = domain.quantities['stage'].vertex_values
855    zv = domain.quantities['elevation'].vertex_values
856    hv = wv - zv
857
858    #Call C-extension
859    from shallow_water_ext import h_limiter_sw
860    hvbar = h_limiter_sw(domain, hc, hv)
861
862    return hvbar
863
864
865def balance_deep_and_shallow(domain):
866    """Compute linear combination between stage as computed by
867    gradient-limiters limiting using w, and stage computed by
868    gradient-limiters limiting using h (h-limiter).
869    The former takes precedence when heights are large compared to the
870    bed slope while the latter takes precedence when heights are
871    relatively small.  Anything in between is computed as a balanced
872    linear combination in order to avoid numerical disturbances which
873    would otherwise appear as a result of hard switching between
874    modes.
875
876    Wrapper for C implementation
877    """
878
879    # FIXME (Ole): I reckon this can be simplified significantly:
880    #
881    # Always use beta_h == 0, and phase it out.
882    # Compute hc and hv in the c-code
883    # Omit updating xmomv
884    #
885    from shallow_water_ext import balance_deep_and_shallow as balance_deep_and_shallow_c
886
887
888    #print 'calling balance depth and shallow'
889    # Shortcuts
890    wc = domain.quantities['stage'].centroid_values
891    zc = domain.quantities['elevation'].centroid_values
892
893    wv = domain.quantities['stage'].vertex_values
894    zv = domain.quantities['elevation'].vertex_values
895
896    # Momentums at centroids
897    xmomc = domain.quantities['xmomentum'].centroid_values
898    ymomc = domain.quantities['ymomentum'].centroid_values
899
900    # Momentums at vertices
901    xmomv = domain.quantities['xmomentum'].vertex_values
902    ymomv = domain.quantities['ymomentum'].vertex_values
903
904    # Limit h
905    if domain.beta_h > 0:
906        hvbar = h_limiter(domain)
907       
908        balance_deep_and_shallow_c(domain, domain.beta_h,
909                                   wc, zc, wv, zv, hvbar,
910                                   xmomc, ymomc, xmomv, ymomv)       
911    else:
912        # print 'Using first order h-limiter'
913        # FIXME: Pass wc in for now - it will be ignored.
914       
915        # This is how one would make a first order h_limited value
916        # as in the old balancer (pre 17 Feb 2005):
917        #  If we wish to hard wire this, one should modify the C-code
918        # from Numeric import zeros, Float
919        # hvbar = zeros( (len(wc), 3), Float)
920        # for i in range(3):
921        #     hvbar[:,i] = wc[:] - zc[:]
922
923        balance_deep_and_shallow_c(domain, domain.beta_h,
924                                   wc, zc, wv, zv, wc, 
925                                   xmomc, ymomc, xmomv, ymomv)
926
927
928
929
930#------------------------------------------------------------------
931# Boundary conditions - specific to the shallow water wave equation
932#------------------------------------------------------------------
933class Reflective_boundary(Boundary):
934    """Reflective boundary returns same conserved quantities as
935    those present in its neighbour volume but reflected.
936
937    This class is specific to the shallow water equation as it
938    works with the momentum quantities assumed to be the second
939    and third conserved quantities.
940    """
941
942    def __init__(self, domain = None):
943        Boundary.__init__(self)
944
945        if domain is None:
946            msg = 'Domain must be specified for reflective boundary'
947            raise msg
948
949        # Handy shorthands
950        self.stage   = domain.quantities['stage'].edge_values
951        self.xmom    = domain.quantities['xmomentum'].edge_values
952        self.ymom    = domain.quantities['ymomentum'].edge_values
953        self.normals = domain.normals
954
955        self.conserved_quantities = zeros(3, Float)
956
957    def __repr__(self):
958        return 'Reflective_boundary'
959
960
961    def evaluate(self, vol_id, edge_id):
962        """Reflective boundaries reverses the outward momentum
963        of the volume they serve.
964        """
965
966        q = self.conserved_quantities
967        q[0] = self.stage[vol_id, edge_id]
968        q[1] = self.xmom[vol_id, edge_id]
969        q[2] = self.ymom[vol_id, edge_id]
970
971        normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
972
973
974        r = rotate(q, normal, direction = 1)
975        r[1] = -r[1]
976        q = rotate(r, normal, direction = -1)
977
978        return q
979
980
981
982class Transmissive_Momentum_Set_Stage_boundary(Boundary):
983    """Returns same momentum conserved quantities as
984    those present in its neighbour volume.
985    Sets stage by specifying a function f of time which may either be a
986    vector function or a scalar function
987
988    Example:
989
990    def waveform(t):
991        return sea_level + normalized_amplitude/cosh(t-25)**2
992
993    Bts = Transmissive_Momentum_Set_Stage_boundary(domain, waveform)
994   
995
996    Underlying domain must be specified when boundary is instantiated
997    """
998
999    def __init__(self, domain = None, function=None):
1000        Boundary.__init__(self)
1001
1002        if domain is None:
1003            msg = 'Domain must be specified for this type boundary'
1004            raise msg
1005
1006        if function is None:
1007            msg = 'Function must be specified for this type boundary'
1008            raise msg
1009
1010        self.domain   = domain
1011        self.function = function
1012
1013    def __repr__(self):
1014        return 'Transmissive_Momentum_Set_Stage_boundary(%s)' %self.domain
1015
1016    def evaluate(self, vol_id, edge_id):
1017        """Transmissive Momentum Set Stage boundaries return the edge momentum
1018        values of the volume they serve.
1019        """
1020
1021        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
1022
1023
1024        t = self.domain.time
1025
1026        if hasattr(self.function, 'time'):
1027            # Roll boundary over if time exceeds           
1028            while t > self.function.time[-1]:
1029                msg = 'WARNING: domain time %.2f has exceeded' %t
1030                msg += 'time provided in '
1031                msg += 'transmissive_momentum_set_stage boundary object.\n'
1032                msg += 'I will continue, reusing the object from t==0'
1033                print msg
1034                t -= self.function.time[-1]
1035
1036
1037        value = self.function(t)
1038        try:
1039            x = float(value)
1040        except:
1041            x = float(value[0])
1042           
1043        q[0] = x
1044        return q
1045
1046
1047        # FIXME: Consider this (taken from File_boundary) to allow
1048        # spatial variation
1049        # if vol_id is not None and edge_id is not None:
1050        #     i = self.boundary_indices[ vol_id, edge_id ]
1051        #     return self.F(t, point_id = i)
1052        # else:
1053        #     return self.F(t)
1054
1055
1056
1057class Dirichlet_Discharge_boundary(Boundary):
1058    """
1059    Sets stage (stage0)
1060    Sets momentum (wh0) in the inward normal direction.
1061
1062    Underlying domain must be specified when boundary is instantiated
1063    """
1064
1065    def __init__(self, domain = None, stage0=None, wh0=None):
1066        Boundary.__init__(self)
1067
1068        if domain is None:
1069            msg = 'Domain must be specified for this type boundary'
1070            raise msg
1071
1072        if stage0 is None:
1073            raise 'set stage'
1074
1075        if wh0 is None:
1076            wh0 = 0.0
1077
1078        self.domain   = domain
1079        self.stage0  = stage0
1080        self.wh0 = wh0
1081
1082    def __repr__(self):
1083        return 'Dirichlet_Discharge_boundary(%s)' %self.domain
1084
1085    def evaluate(self, vol_id, edge_id):
1086        """Set discharge in the (inward) normal direction
1087        """
1088
1089        normal = self.domain.get_normal(vol_id,edge_id)
1090        q = [self.stage0, -self.wh0*normal[0], -self.wh0*normal[1]]
1091        return q
1092
1093
1094        # FIXME: Consider this (taken from File_boundary) to allow
1095        # spatial variation
1096        # if vol_id is not None and edge_id is not None:
1097        #     i = self.boundary_indices[ vol_id, edge_id ]
1098        #     return self.F(t, point_id = i)
1099        # else:
1100        #     return self.F(t)
1101
1102
1103class Field_boundary(Boundary):
1104    """Set boundary from given field represented in an sww file containing values
1105    for stage, xmomentum and ymomentum.
1106    Optionally, the user can specify mean_stage to offset the stage provided in the
1107    sww file.
1108
1109    This function is a thin wrapper around the generic File_boundary. The
1110    difference between the file_boundary and field_boundary is only that the
1111    field_boundary will allow you to change the level of the stage height when
1112    you read in the boundary condition. This is very useful when running
1113    different tide heights in the same area as you need only to convert one
1114    boundary condition to a SWW file, ideally for tide height of 0 m
1115    (saving disk space). Then you can use field_boundary to read this SWW file
1116    and change the stage height (tide) on the fly depending on the scenario.
1117   
1118    """
1119
1120
1121    def __init__(self, filename, domain,
1122                 mean_stage=0.0,
1123                 time_thinning=1, 
1124                 use_cache=False,
1125                 verbose=False):
1126        """Constructor
1127
1128        filename: Name of sww file
1129        domain: pointer to shallow water domain for which the boundary applies
1130        mean_stage: The mean water level which will be added to stage derived
1131                    from the sww file
1132        time_thinning: Will set how many time steps from the sww file read in
1133                       will be interpolated to the boundary. For example if
1134                       the sww file has 1 second time steps and is 24 hours
1135                       in length it has 86400 time steps. If you set
1136                       time_thinning to 1 it will read all these steps.
1137                       If you set it to 100 it will read every 100th step eg
1138                       only 864 step. This parameter is very useful to increase
1139                       the speed of a model run that you are setting up
1140                       and testing.
1141        use_cache:
1142        verbose:
1143       
1144        """
1145
1146        # Create generic file_boundary object
1147        self.file_boundary = File_boundary(filename, domain,
1148                                           time_thinning=time_thinning,
1149                                           use_cache=use_cache,
1150                                           verbose=verbose)
1151       
1152        # Record information from File_boundary
1153        self.F = self.file_boundary.F
1154        self.domain = self.file_boundary.domain
1155       
1156        # Record mean stage
1157        self.mean_stage = mean_stage
1158
1159
1160    def __repr__(self):
1161        return 'Field boundary'
1162
1163
1164    def evaluate(self, vol_id=None, edge_id=None):
1165        """Return linearly interpolated values based on domain.time
1166
1167        vol_id and edge_id are ignored
1168        """
1169
1170        # Evaluate file boundary
1171        q = self.file_boundary.evaluate(vol_id, edge_id)
1172
1173        # Adjust stage
1174        for j, name in enumerate(self.domain.conserved_quantities):
1175            if name == 'stage':
1176                q[j] += self.mean_stage
1177        return q
1178
1179   
1180
1181#-----------------------
1182# Standard forcing terms
1183#-----------------------
1184
1185def gravity(domain):
1186    """Apply gravitational pull in the presence of bed slope
1187    Wrapper calls underlying C implementation
1188    """
1189
1190    xmom = domain.quantities['xmomentum'].explicit_update
1191    ymom = domain.quantities['ymomentum'].explicit_update
1192
1193    stage = domain.quantities['stage']
1194    elevation = domain.quantities['elevation']
1195
1196    h = stage.centroid_values - elevation.centroid_values
1197    z = elevation.vertex_values
1198
1199    x = domain.get_vertex_coordinates()
1200    g = domain.g
1201   
1202
1203    from shallow_water_ext import gravity as gravity_c
1204    gravity_c(g, h, z, x, xmom, ymom) #, 1.0e-6)
1205
1206
1207
1208def manning_friction_implicit(domain):
1209    """Apply (Manning) friction to water momentum   
1210    Wrapper for c version
1211    """
1212
1213
1214    #print 'Implicit friction'
1215
1216    xmom = domain.quantities['xmomentum']
1217    ymom = domain.quantities['ymomentum']
1218
1219    w = domain.quantities['stage'].centroid_values
1220    z = domain.quantities['elevation'].centroid_values
1221
1222    uh = xmom.centroid_values
1223    vh = ymom.centroid_values
1224    eta = domain.quantities['friction'].centroid_values
1225
1226    xmom_update = xmom.semi_implicit_update
1227    ymom_update = ymom.semi_implicit_update
1228
1229    N = len(domain)
1230    eps = domain.minimum_allowed_height
1231    g = domain.g
1232
1233    from shallow_water_ext import manning_friction as manning_friction_c
1234    manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
1235
1236
1237def manning_friction_explicit(domain):
1238    """Apply (Manning) friction to water momentum   
1239    Wrapper for c version
1240    """
1241
1242    # print 'Explicit friction'
1243
1244    xmom = domain.quantities['xmomentum']
1245    ymom = domain.quantities['ymomentum']
1246
1247    w = domain.quantities['stage'].centroid_values
1248    z = domain.quantities['elevation'].centroid_values
1249
1250    uh = xmom.centroid_values
1251    vh = ymom.centroid_values
1252    eta = domain.quantities['friction'].centroid_values
1253
1254    xmom_update = xmom.explicit_update
1255    ymom_update = ymom.explicit_update
1256
1257    N = len(domain)
1258    eps = domain.minimum_allowed_height
1259    g = domain.g
1260
1261    from shallow_water_ext import manning_friction as manning_friction_c
1262    manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
1263
1264
1265# FIXME (Ole): This was implemented for use with one of the analytical solutions (Sampson?)
1266# Is it still needed (30 Oct 2007)?
1267def linear_friction(domain):
1268    """Apply linear friction to water momentum
1269
1270    Assumes quantity: 'linear_friction' to be present
1271    """
1272
1273    from math import sqrt
1274
1275    w = domain.quantities['stage'].centroid_values
1276    z = domain.quantities['elevation'].centroid_values
1277    h = w-z
1278
1279    uh = domain.quantities['xmomentum'].centroid_values
1280    vh = domain.quantities['ymomentum'].centroid_values
1281    tau = domain.quantities['linear_friction'].centroid_values
1282
1283    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1284    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1285
1286    N = len(domain) # number_of_triangles
1287    eps = domain.minimum_allowed_height
1288    g = domain.g #Not necessary? Why was this added?
1289
1290    for k in range(N):
1291        if tau[k] >= eps:
1292            if h[k] >= eps:
1293                S = -tau[k]/h[k]
1294
1295                #Update momentum
1296                xmom_update[k] += S*uh[k]
1297                ymom_update[k] += S*vh[k]
1298
1299
1300
1301#---------------------------------
1302# Experimental auxiliary functions
1303#---------------------------------
1304def check_forcefield(f):
1305    """Check that f is either
1306    1: a callable object f(t,x,y), where x and y are vectors
1307       and that it returns an array or a list of same length
1308       as x and y
1309    2: a scalar
1310    """
1311
1312    if callable(f):
1313        N = 3
1314        x = ones(3, Float)
1315        y = ones(3, Float)
1316        try:
1317            q = f(1.0, x=x, y=y)
1318        except Exception, e:
1319            msg = 'Function %s could not be executed:\n%s' %(f, e)
1320            # FIXME: Reconsider this semantics
1321            raise msg
1322
1323        try:
1324            q = array(q).astype(Float)
1325        except:
1326            msg = 'Return value from vector function %s could ' %f
1327            msg += 'not be converted into a Numeric array of floats.\n'
1328            msg += 'Specified function should return either list or array.'
1329            raise msg
1330
1331        # Is this really what we want?
1332        msg = 'Return vector from function %s ' %f
1333        msg += 'must have same lenght as input vectors'
1334        assert len(q) == N, msg
1335
1336    else:
1337        try:
1338            f = float(f)
1339        except:
1340            msg = 'Force field %s must be either a scalar' %f
1341            msg += ' or a vector function'
1342            raise Exception(msg)
1343    return f
1344
1345
1346class Wind_stress:
1347    """Apply wind stress to water momentum in terms of
1348    wind speed [m/s] and wind direction [degrees]
1349    """
1350
1351    def __init__(self, *args, **kwargs):
1352        """Initialise windfield from wind speed s [m/s]
1353        and wind direction phi [degrees]
1354
1355        Inputs v and phi can be either scalars or Python functions, e.g.
1356
1357        W = Wind_stress(10, 178)
1358
1359        #FIXME - 'normal' degrees are assumed for now, i.e. the
1360        vector (1,0) has zero degrees.
1361        We may need to convert from 'compass' degrees later on and also
1362        map from True north to grid north.
1363
1364        Arguments can also be Python functions of t,x,y as in
1365
1366        def speed(t,x,y):
1367            ...
1368            return s
1369
1370        def angle(t,x,y):
1371            ...
1372            return phi
1373
1374        where x and y are vectors.
1375
1376        and then pass the functions in
1377
1378        W = Wind_stress(speed, angle)
1379
1380        The instantiated object W can be appended to the list of
1381        forcing_terms as in
1382
1383        Alternatively, one vector valued function for (speed, angle)
1384        can be applied, providing both quantities simultaneously.
1385        As in
1386        W = Wind_stress(F), where returns (speed, angle) for each t.
1387
1388        domain.forcing_terms.append(W)
1389        """
1390
1391        from anuga.config import rho_a, rho_w, eta_w
1392        from Numeric import array, Float
1393
1394        if len(args) == 2:
1395            s = args[0]
1396            phi = args[1]
1397        elif len(args) == 1:
1398            # Assume vector function returning (s, phi)(t,x,y)
1399            vector_function = args[0]
1400            s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
1401            phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
1402        else:
1403           # Assume info is in 2 keyword arguments
1404
1405           if len(kwargs) == 2:
1406               s = kwargs['s']
1407               phi = kwargs['phi']
1408           else:
1409               raise 'Assumes two keyword arguments: s=..., phi=....'
1410
1411        self.speed = check_forcefield(s)
1412        self.phi = check_forcefield(phi)
1413
1414        self.const = eta_w*rho_a/rho_w
1415
1416
1417    def __call__(self, domain):
1418        """Evaluate windfield based on values found in domain
1419        """
1420
1421        from math import pi, cos, sin, sqrt
1422        from Numeric import Float, ones, ArrayType
1423
1424        xmom_update = domain.quantities['xmomentum'].explicit_update
1425        ymom_update = domain.quantities['ymomentum'].explicit_update
1426
1427        N = len(domain) # number_of_triangles
1428        t = domain.time
1429
1430        if callable(self.speed):
1431            xc = domain.get_centroid_coordinates()
1432            s_vec = self.speed(t, xc[:,0], xc[:,1])
1433        else:
1434            # Assume s is a scalar
1435
1436            try:
1437                s_vec = self.speed * ones(N, Float)
1438            except:
1439                msg = 'Speed must be either callable or a scalar: %s' %self.s
1440                raise msg
1441
1442
1443        if callable(self.phi):
1444            xc = domain.get_centroid_coordinates()
1445            phi_vec = self.phi(t, xc[:,0], xc[:,1])
1446        else:
1447            # Assume phi is a scalar
1448
1449            try:
1450                phi_vec = self.phi * ones(N, Float)
1451            except:
1452                msg = 'Angle must be either callable or a scalar: %s' %self.phi
1453                raise msg
1454
1455        assign_windfield_values(xmom_update, ymom_update,
1456                                s_vec, phi_vec, self.const)
1457
1458
1459def assign_windfield_values(xmom_update, ymom_update,
1460                            s_vec, phi_vec, const):
1461    """Python version of assigning wind field to update vectors.
1462    A c version also exists (for speed)
1463    """
1464    from math import pi, cos, sin, sqrt
1465
1466    N = len(s_vec)
1467    for k in range(N):
1468        s = s_vec[k]
1469        phi = phi_vec[k]
1470
1471        # Convert to radians
1472        phi = phi*pi/180
1473
1474        # Compute velocity vector (u, v)
1475        u = s*cos(phi)
1476        v = s*sin(phi)
1477
1478        # Compute wind stress
1479        S = const * sqrt(u**2 + v**2)
1480        xmom_update[k] += S*u
1481        ymom_update[k] += S*v
1482
1483
1484
1485
1486
1487class General_forcing:
1488    """Class General_forcing - general explicit forcing term for update of quantity
1489   
1490    This is used by Inflow and Rainfall for instance
1491   
1492
1493    General_forcing(quantity_name, rate, center, radius, polygon)
1494
1495    domain:     ANUGA computational domain
1496    quantity_name: Name of quantity to update. It must be a known conserved quantity.
1497    rate [?/s]: Total rate of change over the specified area. 
1498                This parameter can be either a constant or a
1499                function of time. Positive values indicate increases,
1500                negative values indicate decreases.
1501                Rate can be None at initialisation but must be specified
1502                before forcting term is applied (i.e. simulation has started).
1503
1504    center [m]: Coordinates at center of flow point
1505    radius [m]: Size of circular area
1506    polygon:    Arbitrary polygon.
1507
1508
1509    Either center, radius or polygon can be specified but not both.
1510    If neither are specified the entire domain gets updated.
1511   
1512    See Inflow or Rainfall for examples of use
1513    """
1514
1515
1516    # FIXME (AnyOne) : Add various methods to allow spatial variations
1517
1518    def __init__(self,
1519                 domain,
1520                 quantity_name,
1521                 rate=None,
1522                 center=None, radius=None,
1523                 polygon=None,
1524                 verbose=False):
1525                     
1526
1527        from math import pi
1528
1529        self.domain = domain
1530        self.quantity_name = quantity_name
1531        self.rate = rate
1532        self.center = ensure_numeric(center)
1533        self.radius = radius
1534        self.polygon = polygon       
1535        self.verbose = verbose
1536
1537        # Update area if applicable
1538        self.area = None       
1539        if center is not None and radius is not None:
1540            assert len(center) == 2
1541            msg = 'Polygon cannot be specified when center and radius are'
1542            assert polygon is None, msg
1543
1544            self.area = radius**2*pi
1545       
1546        if polygon is not None:
1547            self.area = polygon_area(self.polygon)
1548
1549
1550        # Pointer to update vector
1551        self.update = domain.quantities[self.quantity_name].explicit_update           
1552
1553        # Determine indices in flow area
1554        N = len(domain)   
1555        points = domain.get_centroid_coordinates(absolute=True)
1556
1557        self.indices = None
1558        if self.center is not None and self.radius is not None:
1559            # Inlet is circular
1560           
1561            self.indices = []
1562            for k in range(N):
1563                x, y = points[k,:] # Centroid
1564                if ((x-self.center[0])**2+(y-self.center[1])**2) < self.radius**2:
1565                    self.indices.append(k)
1566                   
1567        if self.polygon is not None:                   
1568            # Inlet is polygon
1569            self.indices = inside_polygon(points, self.polygon)
1570           
1571
1572           
1573
1574
1575    def __call__(self, domain):
1576        """Apply inflow function at time specified in domain and update stage
1577        """
1578
1579        # Call virtual method allowing local modifications
1580        rate = self.update_rate(domain.get_time())
1581        if rate is None:
1582            msg = 'Attribute rate must be specified in General_forcing'
1583            msg += ' or its descendants before attempting to call it'
1584            raise Exception, msg
1585       
1586
1587        # Now rate is a number
1588        if self.verbose is True:
1589            print 'Rate of %s at time = %.2f = %f' %(self.quantity_name,
1590                                                     domain.get_time(),
1591                                                     rate)
1592
1593
1594        if self.indices is None:
1595            self.update[:] += rate
1596        else:
1597            # Brute force assignment of restricted rate
1598            for k in self.indices:
1599                self.update[k] += rate
1600
1601
1602    def update_rate(self, t):
1603        """Virtual method allowing local modifications by writing an
1604        overriding version in descendant
1605       
1606        """
1607        if callable(self.rate):
1608            rate = self.rate(t)
1609        else:
1610            rate = self.rate
1611
1612        return rate
1613
1614
1615    def get_quantity_values(self):
1616        """Return values for specified quantity restricted to opening
1617        """
1618        return self.domain.quantities[self.quantity_name].get_values(indices=self.indices)
1619   
1620
1621    def set_quantity_values(self, val):
1622        """Set values for specified quantity restricted to opening
1623        """
1624        self.domain.quantities[self.quantity_name].set_values(val, indices=self.indices)   
1625
1626
1627
1628class Rainfall(General_forcing):
1629    """Class Rainfall - general 'rain over entire domain' forcing term.
1630   
1631    Used for implementing Rainfall over the entire domain.
1632       
1633        Current Limited to only One Gauge..
1634       
1635        Need to add Spatial Varying Capability
1636        (This module came from copying and amending the Inflow Code)
1637   
1638    Rainfall(rain)
1639
1640    domain   
1641    rain [mm/s]:  Total rain rate over the specified domain. 
1642                  NOTE: Raingauge Data needs to reflect the time step.
1643                  IE: if Gauge is mm read at a time step, then the input
1644                  here is as mm/(timeStep) so 10mm in 5minutes becomes
1645                  10/(5x60) = 0.0333mm/s.
1646       
1647       
1648                  This parameter can be either a constant or a
1649                  function of time. Positive values indicate inflow,
1650                  negative values indicate outflow.
1651                  (and be used for Infiltration - Write Seperate Module)
1652                  The specified flow will be divided by the area of
1653                  the inflow region and then applied to update the
1654                  stage quantity.
1655
1656    polygon: Specifies a polygon to restrict the rainfall.
1657   
1658    Examples
1659    How to put them in a run File...
1660       
1661    #--------------------------------------------------------------------------
1662    # Setup specialised forcing terms
1663    #--------------------------------------------------------------------------
1664    # This is the new element implemented by Ole and Rudy to allow direct
1665    # input of Inflow in mm/s
1666
1667    catchmentrainfall = Rainfall(rain=file_function('Q100_2hr_Rain.tms')) 
1668                        # Note need path to File in String.
1669                        # Else assumed in same directory
1670
1671    domain.forcing_terms.append(catchmentrainfall)
1672    """
1673
1674   
1675    def __init__(self,
1676                 domain,
1677                 rate=0.0,
1678                 center=None, radius=None,
1679                 polygon=None,
1680                 verbose=False):
1681
1682        # Converting mm/s to m/s to apply in ANUGA)
1683        if callable(rate):
1684            rain = lambda t: rate(t)/1000.0
1685        else:
1686            rain = rate/1000.0           
1687           
1688        General_forcing.__init__(self,
1689                                 domain,
1690                                 'stage',
1691                                 rate=rain,
1692                                 center=center, radius=radius,
1693                                 polygon=polygon,
1694                                 verbose=verbose)
1695
1696       
1697
1698
1699
1700
1701class Inflow(General_forcing):
1702    """Class Inflow - general 'rain and drain' forcing term.
1703   
1704    Useful for implementing flows in and out of the domain.
1705   
1706    Inflow(flow, center, radius, polygon)
1707
1708    domain
1709    flow [m^3/s]: Total flow rate over the specified area. 
1710                  This parameter can be either a constant or a
1711                  function of time. Positive values indicate inflow,
1712                  negative values indicate outflow.
1713                  The specified flow will be divided by the area of
1714                  the inflow region and then applied to update the
1715                  quantity in question.     
1716    center [m]: Coordinates at center of flow point
1717    radius [m]: Size of circular area
1718    polygon:    Arbitrary polygon.
1719
1720    Either center, radius or polygon must be specified
1721   
1722    Examples
1723
1724    # Constant drain at 0.003 m^3/s.
1725    # The outflow area is 0.07**2*pi=0.0154 m^2
1726    # This corresponds to a rate of change of 0.003/0.0154 = 0.2 m/s
1727    #                                     
1728    Inflow((0.7, 0.4), 0.07, -0.003)
1729
1730
1731    # Tap turning up to a maximum inflow of 0.0142 m^3/s.
1732    # The inflow area is 0.03**2*pi = 0.00283 m^2
1733    # This corresponds to a rate of change of 0.0142/0.00283 = 5 m/s     
1734    # over the specified area
1735    Inflow((0.5, 0.5), 0.03, lambda t: min(0.01*t, 0.0142))
1736
1737
1738    #--------------------------------------------------------------------------
1739    # Setup specialised forcing terms
1740    #--------------------------------------------------------------------------
1741    # This is the new element implemented by Ole to allow direct input
1742    # of Inflow in m^3/s
1743
1744    hydrograph = Inflow(center=(320, 300), radius=10,
1745                        flow=file_function('Q/QPMF_Rot_Sub13.tms'))
1746
1747    domain.forcing_terms.append(hydrograph)
1748   
1749    """
1750
1751
1752    def __init__(self,
1753                 domain,
1754                 rate=0.0,
1755                 center=None, radius=None,
1756                 polygon=None,
1757                 verbose=False):                 
1758
1759
1760        #msg = 'Class Inflow must have either center & radius or a polygon specified.'
1761        #assert center is not None and radius is not None or\
1762        #       polygon is not None, msg
1763
1764
1765        # Create object first to make area is available
1766        General_forcing.__init__(self,
1767                                 domain,
1768                                 'stage',
1769                                 rate=rate,
1770                                 center=center, radius=radius,
1771                                 polygon=polygon,
1772                                 verbose=verbose)
1773
1774    def update_rate(self, t):
1775        """Virtual method allowing local modifications by writing an
1776        overriding version in descendant
1777
1778        This one converts m^3/s to m/s which can be added directly to 'stage' in ANUGA
1779        """
1780
1781       
1782       
1783        if callable(self.rate):
1784            _rate = self.rate(t)/self.area
1785        else:
1786            _rate = self.rate/self.area
1787
1788        return _rate
1789
1790
1791
1792
1793#------------------
1794# Initialise module
1795#------------------
1796
1797
1798from anuga.utilities import compile
1799if compile.can_use_C_extension('shallow_water_ext.c'):
1800    # Underlying C implementations can be accessed
1801
1802    from shallow_water_ext import rotate, assign_windfield_values
1803else:
1804    msg = 'C implementations could not be accessed by %s.\n ' %__file__
1805    msg += 'Make sure compile_all.py has been run as described in '
1806    msg += 'the ANUGA installation guide.'
1807    raise Exception, msg
1808
1809
1810# Optimisation with psyco
1811from anuga.config import use_psyco
1812if use_psyco:
1813    try:
1814        import psyco
1815    except:
1816        import os
1817        if os.name == 'posix' and os.uname()[4] == 'x86_64':
1818            pass
1819            #Psyco isn't supported on 64 bit systems, but it doesn't matter
1820        else:
1821            msg = 'WARNING: psyco (speedup) could not import'+\
1822                  ', you may want to consider installing it'
1823            print msg
1824    else:
1825        psyco.bind(Domain.distribute_to_vertices_and_edges)
1826        psyco.bind(Domain.compute_fluxes)
1827
1828if __name__ == "__main__":
1829    pass
1830
1831
Note: See TracBrowser for help on using the repository browser.