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

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

Added input tests regarding polygons and points

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