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

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

Implemented default_boundary option in File_boundary and Field_boundary as
per ticket:293 and added a note in the user manual.

  • Property svn:keywords set to LastChangedDate LastChangedRevision LastChangedBy HeadURL Id
File size: 60.3 KB
RevLine 
[4004]1"""Finite-volume computations of the shallow water wave equation.
2
[4005]3Title: ANGUA shallow_water_domain - 2D triangular domains for finite-volume
4       computations of the shallow water wave equation.
[3804]5
6
[5186]7Author: Ole Nielsen, Ole.Nielsen@ga.gov.au
8        Stephen Roberts, Stephen.Roberts@anu.edu.au
9        Duncan Gray, Duncan.Gray@ga.gov.au
[4004]10
11CreationDate: 2004
12
13Description:
[4005]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
[4004]17
[4005]18    U_t + E_x + G_y = S
[3804]19
[4005]20    where
[3804]21
[4005]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, ...)
[3804]27
[4005]28    and _t, _x, _y denote the derivative with respect to t, x and y
29    respectively.
[3804]30
31
[4005]32    The quantities are
[3804]33
[4005]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]
[3804]44
[4005]45    eta                        mannings friction coefficient [to appear]
46    nu                         wind stress coefficient [to appear]
[3804]47
[4005]48    The conserved quantities are w, uh, vh
[3804]49
[4004]50Reference:
[4005]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
[3804]54
[4005]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
[3804]61
62
[4005]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 $
[3804]67
[4004]68Constraints: See GPL license in the user guide
69Version: 1.0 ($Revision: 5657 $)
70ModifiedBy:
[4005]71    $Author: ole $
72    $Date: 2008-08-14 00:26:06 +0000 (Thu, 14 Aug 2008) $
[4004]73
[3804]74"""
75
[4769]76# Subversion keywords:
[3804]77#
[4769]78# $LastChangedDate: 2008-08-14 00:26:06 +0000 (Thu, 14 Aug 2008) $
79# $LastChangedRevision: 5657 $
80# $LastChangedBy: ole $
[3804]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
[5294]98from anuga.utilities.numerical_tools import gradient, mean, ensure_numeric
[3804]99from anuga.config import minimum_storable_height
100from anuga.config import minimum_allowed_height, maximum_allowed_speed
[5442]101from anuga.config import g, epsilon, beta_w, beta_w_dry,\
[4631]102     beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, tight_slope_limiters
[3876]103from anuga.config import alpha_balance
[4685]104from anuga.config import optimise_dry_cells
[5162]105from anuga.config import optimised_gradient_limiter
[5175]106from anuga.config import use_edge_limiter
[5290]107from anuga.config import use_centroid_velocities
[3804]108
[5290]109
[5570]110from anuga.utilities.polygon import inside_polygon, polygon_area, is_inside_polygon
[5294]111
112
113
[4769]114#---------------------
115# Shallow water domain
116#---------------------
[3804]117class Domain(Generic_Domain):
118
[4454]119    conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
120    other_quantities = ['elevation', 'friction']
121   
[3804]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,
[3926]135                 numproc=1,
[3928]136                 number_of_full_nodes=None,
137                 number_of_full_triangles=None):
[3804]138
139
140        other_quantities = ['elevation', 'friction']
141        Generic_Domain.__init__(self,
142                                coordinates,
143                                vertices,
144                                boundary,
[4454]145                                Domain.conserved_quantities,
146                                Domain.other_quantities,
[3804]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,
[3926]156                                numproc,
157                                number_of_full_nodes=number_of_full_nodes,
158                                number_of_full_triangles=number_of_full_triangles) 
[3804]159
[4769]160
[4258]161        self.set_minimum_allowed_height(minimum_allowed_height)
162       
[3804]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
[3876]171        self.alpha_balance = alpha_balance
[3804]172
[4631]173        self.tight_slope_limiters = tight_slope_limiters
[4685]174        self.optimise_dry_cells = optimise_dry_cells
[4239]175
[4769]176        self.forcing_terms.append(manning_friction_implicit)
[3804]177        self.forcing_terms.append(gravity)
178
[4701]179        # Stored output
[3804]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']
[5162]185
186        # Limiters
[5175]187        self.use_edge_limiter = use_edge_limiter
[5162]188        self.optimised_gradient_limiter = optimised_gradient_limiter
[5290]189        self.use_centroid_velocities = use_centroid_velocities
[3804]190
[5290]191
[3848]192    def set_all_limiters(self, beta):
[3847]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        """
[3804]197
[3847]198        self.beta_w      = beta
199        self.beta_w_dry  = beta
[5162]200        self.quantities['stage'].beta = beta
201       
[3847]202        self.beta_uh     = beta
203        self.beta_uh_dry = beta
[5162]204        self.quantities['xmomentum'].beta = beta
205       
[3847]206        self.beta_vh     = beta
207        self.beta_vh_dry = beta
[5162]208        self.quantities['ymomentum'].beta = beta
209       
[3847]210       
211
[3804]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        """
[3954]217
218        # FIXME (Ole): how about using the word continuous vertex values?
[3804]219        self.smooth = not flag
220
[4733]221        # Reduction operation for get_vertex_values
[3804]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
[4258]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        """
[4438]248
249        #FIXME (Ole): rename H0 to minimum_allowed_height_in_flux_computation
[4701]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.
[4258]254        self.minimum_allowed_height = minimum_allowed_height
255        self.H0 = minimum_allowed_height   
256       
[3804]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
[4254]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       
[4701]279       
[3804]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
[4701]304        # Check correcness
[3804]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):
[4554]362        """Return location of highest elevation where h > 0
[3804]363
364        Optional argument:
365            indices is the set of element ids that the operation applies to.
366
367        Usage:
[4554]368            q = get_maximum_inundation_location()
[3804]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):
[4769]400        # Call correct module function
[5176]401        if self.use_edge_limiter:
[5306]402            distribute_using_edge_limiter(self)           
[5175]403        else:
[5306]404            distribute_using_vertex_limiter(self)
[3804]405
406
[5175]407
408
[3804]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
[4769]417        # Call check integrity here rather than from user scripts
418        # self.check_integrity()
[3804]419
[5162]420        msg = 'Parameter beta_w must be in the interval [0, 2['
421        assert 0 <= self.beta_w <= 2.0, msg
[3804]422
423
[4769]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
[3804]428        self.distribute_to_vertices_and_edges()
429
430        if self.store is True and self.time == 0.0:
431            self.initialise_storage()
[4769]432            # print 'Storing results in ' + self.writer.filename
[3804]433        else:
434            pass
[4769]435            # print 'Results will not be stored.'
436            # print 'To store results set domain.store = True'
437            # FIXME: Diagnostic output should be controlled by
[3804]438            # a 'verbose' flag living in domain (or in a parent class)
439
[4769]440        # Call basic machinery from parent class
[3804]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
[4769]447            # Store model data, e.g. for subsequent visualisation
[3804]448            if self.store is True:
449                self.store_timestep(self.quantities_to_be_stored)
450
[4769]451            # FIXME: Could maybe be taken from specified list
452            # of 'store every step' quantities
[3804]453
[4769]454            # Pass control on to outer loop for more specific actions
[3804]455            yield(t)
456
[4769]457
[3804]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
[4769]465        # Initialise writer
[3804]466        self.writer = get_dataobject(self, mode = 'w')
467
[4769]468        # Store vertices and connectivity
[3804]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
[4827]480       
[4836]481    def timestepping_statistics(self,
482                                track_speeds=False,
483                                triangle_id=None):       
[4827]484        """Return string with time stepping statistics for printing or logging
[3804]485
[4827]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
[4836]496        msg = Generic_Domain.timestepping_statistics(self,
497                                                     track_speeds,
498                                                     triangle_id)
[4827]499
500        if track_speeds is True:
501
502            # qwidth determines the text field used for quantities
503            qwidth = self.qwidth
504       
[4836]505            # Selected triangle
[4827]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
[4769]622#=============== End of class Shallow Water Domain ===============================
[3804]623
624
[4769]625#-----------------
[3804]626# Flux computation
[4769]627#-----------------
[3804]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.
[4769]647   
648
649    This wrapper calls the underlying C version of compute fluxes
[3804]650    """
651
652    import sys
653
[3928]654    N = len(domain) # number_of_triangles
[3804]655
[4769]656    # Shortcuts
[3804]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)
[4769]663    from shallow_water_ext import\
664         compute_fluxes_ext_central as compute_fluxes_ext
[3804]665
666
[4769]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))
[3804]691
[4769]692    domain.flux_timestep = flux_timestep
[3804]693
694
695
[4769]696#---------------------------------------
697# Module functions for gradient limiting
698#---------------------------------------
[3804]699
700
[4769]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):
[3804]704    """Wrapper calling C version of extrapolate_second_order_sw
705    """
706    import sys
707
[3928]708    N = len(domain) # number_of_triangles
[3804]709
[4710]710    # Shortcuts
[3804]711    Stage = domain.quantities['stage']
712    Xmom = domain.quantities['xmomentum']
713    Ymom = domain.quantities['ymomentum']
714    Elevation = domain.quantities['elevation']
[4710]715
[4769]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,
[5315]730              int(domain.optimise_dry_cells))
[3804]731
732
[5306]733def distribute_using_vertex_limiter(domain):
[3804]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
[5162]756   
[3804]757
[4769]758    # Remove very thin layers of water
[3804]759    protect_against_infinitesimal_and_negative_heights(domain)
760
[4769]761    # Extrapolate all conserved quantities
[5162]762    if domain.optimised_gradient_limiter:
[4769]763        # MH090605 if second order,
764        # perform the extrapolation and limiting on
765        # all of the conserved quantities
[3804]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:
[4769]776        # Old code:
[3804]777        for name in domain.conserved_quantities:
778            Q = domain.quantities[name]
[4701]779
[3804]780            if domain._order_ == 1:
781                Q.extrapolate_first_order()
782            elif domain._order_ == 2:
[5306]783                Q.extrapolate_second_order_and_limit_by_vertex()
[3804]784            else:
785                raise 'Unknown order'
786
787
[5290]788    # Take bed elevation into account when water heights are small
[3804]789    balance_deep_and_shallow(domain)
790
[5290]791    # Compute edge values by interpolation
[3804]792    for name in domain.conserved_quantities:
793        Q = domain.quantities[name]
794        Q.interpolate_from_vertices_to_edges()
795
796
[5306]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
[3804]841def protect_against_infinitesimal_and_negative_heights(domain):
842    """Protect against infinitesimal heights and associated high velocities
843    """
844
[4769]845    # Shortcuts
[3804]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
[4769]869    Wrapper for C implementation
[3804]870    """
871
[4769]872    from shallow_water_ext import balance_deep_and_shallow as balance_deep_and_shallow_c
[5175]873
874
[4733]875    # Shortcuts
[3804]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
[4733]882    # Momentums at centroids
[3804]883    xmomc = domain.quantities['xmomentum'].centroid_values
884    ymomc = domain.quantities['ymomentum'].centroid_values
885
[4733]886    # Momentums at vertices
[3804]887    xmomv = domain.quantities['xmomentum'].vertex_values
888    ymomv = domain.quantities['ymomentum'].vertex_values
889
[5442]890    balance_deep_and_shallow_c(domain,
891                               wc, zc, wv, zv, wc, 
892                               xmomc, ymomc, xmomv, ymomv)
[3804]893
894
895
896
[4769]897#------------------------------------------------------------------
898# Boundary conditions - specific to the shallow water wave equation
899#------------------------------------------------------------------
[3804]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
[4769]916        # Handy shorthands
[3804]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
[5089]990
[5081]991        t = self.domain.time
992
[5089]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]
[5081]1002
[5089]1003
[5081]1004        value = self.function(t)
[3804]1005        try:
1006            x = float(value)
[5081]1007        except:
[3804]1008            x = float(value[0])
1009           
1010        q[0] = x
1011        return q
1012
1013
[4769]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)
[3804]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
[4769]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)
[3804]1068
1069
[4312]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.
[3804]1075
[4754]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   
[4312]1085    """
[3804]1086
1087
[4312]1088    def __init__(self, filename, domain,
1089                 mean_stage=0.0,
[5657]1090                 time_thinning=1,
1091                 boundary_polygon=None,   
1092                 default_boundary=None,                 
[4312]1093                 use_cache=False,
[5657]1094                 verbose=False):
[4312]1095        """Constructor
1096
1097        filename: Name of sww file
1098        domain: pointer to shallow water domain for which the boundary applies
[4754]1099        mean_stage: The mean water level which will be added to stage derived
1100                    from the sww file
1101        time_thinning: Will set how many time steps from the sww file read in
1102                       will be interpolated to the boundary. For example if
1103                       the sww file has 1 second time steps and is 24 hours
1104                       in length it has 86400 time steps. If you set
1105                       time_thinning to 1 it will read all these steps.
1106                       If you set it to 100 it will read every 100th step eg
1107                       only 864 step. This parameter is very useful to increase
1108                       the speed of a model run that you are setting up
1109                       and testing.
[5657]1110                       
1111        default_boundary: Must be either None or an instance of a
1112                          class descending from class Boundary.
1113                          This will be used in case model time exceeds
1114                          that available in the underlying data.
1115                                               
[4312]1116        use_cache:
1117        verbose:
1118       
1119        """
1120
1121        # Create generic file_boundary object
1122        self.file_boundary = File_boundary(filename, domain,
1123                                           time_thinning=time_thinning,
[5657]1124                                           boundary_polygon=boundary_polygon,
1125                                           default_boundary=default_boundary,
[4312]1126                                           use_cache=use_cache,
[5657]1127                                           verbose=verbose)
1128
[4312]1129       
1130        # Record information from File_boundary
1131        self.F = self.file_boundary.F
1132        self.domain = self.file_boundary.domain
1133       
1134        # Record mean stage
1135        self.mean_stage = mean_stage
1136
1137
1138    def __repr__(self):
1139        return 'Field boundary'
1140
1141
1142    def evaluate(self, vol_id=None, edge_id=None):
1143        """Return linearly interpolated values based on domain.time
1144
1145        vol_id and edge_id are ignored
1146        """
1147
1148        # Evaluate file boundary
1149        q = self.file_boundary.evaluate(vol_id, edge_id)
1150
1151        # Adjust stage
1152        for j, name in enumerate(self.domain.conserved_quantities):
1153            if name == 'stage':
1154                q[j] += self.mean_stage
1155        return q
1156
1157   
1158
[4769]1159#-----------------------
1160# Standard forcing terms
1161#-----------------------
1162
[3804]1163def gravity(domain):
1164    """Apply gravitational pull in the presence of bed slope
[4769]1165    Wrapper calls underlying C implementation
[3804]1166    """
1167
1168    xmom = domain.quantities['xmomentum'].explicit_update
1169    ymom = domain.quantities['ymomentum'].explicit_update
1170
[4687]1171    stage = domain.quantities['stage']
1172    elevation = domain.quantities['elevation']
[3804]1173
[4687]1174    h = stage.centroid_values - elevation.centroid_values
1175    z = elevation.vertex_values
1176
[3804]1177    x = domain.get_vertex_coordinates()
1178    g = domain.g
[4687]1179   
[3804]1180
[4769]1181    from shallow_water_ext import gravity as gravity_c
1182    gravity_c(g, h, z, x, xmom, ymom) #, 1.0e-6)
[3804]1183
1184
[4438]1185
[4769]1186def manning_friction_implicit(domain):
1187    """Apply (Manning) friction to water momentum   
1188    Wrapper for c version
[3804]1189    """
1190
1191
1192    #print 'Implicit friction'
1193
1194    xmom = domain.quantities['xmomentum']
1195    ymom = domain.quantities['ymomentum']
1196
1197    w = domain.quantities['stage'].centroid_values
1198    z = domain.quantities['elevation'].centroid_values
1199
1200    uh = xmom.centroid_values
1201    vh = ymom.centroid_values
1202    eta = domain.quantities['friction'].centroid_values
1203
1204    xmom_update = xmom.semi_implicit_update
1205    ymom_update = ymom.semi_implicit_update
1206
[3928]1207    N = len(domain)
[3804]1208    eps = domain.minimum_allowed_height
1209    g = domain.g
1210
[4769]1211    from shallow_water_ext import manning_friction as manning_friction_c
1212    manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
[3804]1213
1214
[4769]1215def manning_friction_explicit(domain):
1216    """Apply (Manning) friction to water momentum   
1217    Wrapper for c version
[3804]1218    """
1219
[4769]1220    # print 'Explicit friction'
[3804]1221
1222    xmom = domain.quantities['xmomentum']
1223    ymom = domain.quantities['ymomentum']
1224
1225    w = domain.quantities['stage'].centroid_values
1226    z = domain.quantities['elevation'].centroid_values
1227
1228    uh = xmom.centroid_values
1229    vh = ymom.centroid_values
1230    eta = domain.quantities['friction'].centroid_values
1231
1232    xmom_update = xmom.explicit_update
1233    ymom_update = ymom.explicit_update
1234
[3928]1235    N = len(domain)
[3804]1236    eps = domain.minimum_allowed_height
1237    g = domain.g
1238
[4769]1239    from shallow_water_ext import manning_friction as manning_friction_c
1240    manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
[3804]1241
1242
[4769]1243# FIXME (Ole): This was implemented for use with one of the analytical solutions (Sampson?)
1244# Is it still needed (30 Oct 2007)?
[3804]1245def linear_friction(domain):
1246    """Apply linear friction to water momentum
1247
1248    Assumes quantity: 'linear_friction' to be present
1249    """
1250
1251    from math import sqrt
1252
1253    w = domain.quantities['stage'].centroid_values
1254    z = domain.quantities['elevation'].centroid_values
1255    h = w-z
1256
1257    uh = domain.quantities['xmomentum'].centroid_values
1258    vh = domain.quantities['ymomentum'].centroid_values
1259    tau = domain.quantities['linear_friction'].centroid_values
1260
1261    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1262    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1263
[3928]1264    N = len(domain) # number_of_triangles
[3804]1265    eps = domain.minimum_allowed_height
1266    g = domain.g #Not necessary? Why was this added?
1267
1268    for k in range(N):
1269        if tau[k] >= eps:
1270            if h[k] >= eps:
1271                S = -tau[k]/h[k]
1272
1273                #Update momentum
1274                xmom_update[k] += S*uh[k]
1275                ymom_update[k] += S*vh[k]
1276
1277
1278
[4769]1279#---------------------------------
1280# Experimental auxiliary functions
1281#---------------------------------
[3804]1282def check_forcefield(f):
1283    """Check that f is either
1284    1: a callable object f(t,x,y), where x and y are vectors
1285       and that it returns an array or a list of same length
1286       as x and y
1287    2: a scalar
1288    """
1289
1290    if callable(f):
1291        N = 3
1292        x = ones(3, Float)
1293        y = ones(3, Float)
1294        try:
1295            q = f(1.0, x=x, y=y)
1296        except Exception, e:
1297            msg = 'Function %s could not be executed:\n%s' %(f, e)
[4769]1298            # FIXME: Reconsider this semantics
[3804]1299            raise msg
1300
1301        try:
1302            q = array(q).astype(Float)
1303        except:
1304            msg = 'Return value from vector function %s could ' %f
1305            msg += 'not be converted into a Numeric array of floats.\n'
1306            msg += 'Specified function should return either list or array.'
1307            raise msg
1308
[4769]1309        # Is this really what we want?
[3804]1310        msg = 'Return vector from function %s ' %f
1311        msg += 'must have same lenght as input vectors'
1312        assert len(q) == N, msg
1313
1314    else:
1315        try:
1316            f = float(f)
1317        except:
1318            msg = 'Force field %s must be either a scalar' %f
1319            msg += ' or a vector function'
1320            raise Exception(msg)
1321    return f
1322
1323
1324class Wind_stress:
1325    """Apply wind stress to water momentum in terms of
1326    wind speed [m/s] and wind direction [degrees]
1327    """
1328
1329    def __init__(self, *args, **kwargs):
1330        """Initialise windfield from wind speed s [m/s]
1331        and wind direction phi [degrees]
1332
1333        Inputs v and phi can be either scalars or Python functions, e.g.
1334
1335        W = Wind_stress(10, 178)
1336
1337        #FIXME - 'normal' degrees are assumed for now, i.e. the
1338        vector (1,0) has zero degrees.
1339        We may need to convert from 'compass' degrees later on and also
1340        map from True north to grid north.
1341
1342        Arguments can also be Python functions of t,x,y as in
1343
1344        def speed(t,x,y):
1345            ...
1346            return s
1347
1348        def angle(t,x,y):
1349            ...
1350            return phi
1351
1352        where x and y are vectors.
1353
1354        and then pass the functions in
1355
1356        W = Wind_stress(speed, angle)
1357
1358        The instantiated object W can be appended to the list of
1359        forcing_terms as in
1360
1361        Alternatively, one vector valued function for (speed, angle)
1362        can be applied, providing both quantities simultaneously.
1363        As in
1364        W = Wind_stress(F), where returns (speed, angle) for each t.
1365
1366        domain.forcing_terms.append(W)
1367        """
1368
1369        from anuga.config import rho_a, rho_w, eta_w
1370        from Numeric import array, Float
1371
1372        if len(args) == 2:
1373            s = args[0]
1374            phi = args[1]
1375        elif len(args) == 1:
[4769]1376            # Assume vector function returning (s, phi)(t,x,y)
[3804]1377            vector_function = args[0]
1378            s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
1379            phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
1380        else:
[4769]1381           # Assume info is in 2 keyword arguments
[3804]1382
1383           if len(kwargs) == 2:
1384               s = kwargs['s']
1385               phi = kwargs['phi']
1386           else:
1387               raise 'Assumes two keyword arguments: s=..., phi=....'
1388
1389        self.speed = check_forcefield(s)
1390        self.phi = check_forcefield(phi)
1391
1392        self.const = eta_w*rho_a/rho_w
1393
1394
1395    def __call__(self, domain):
1396        """Evaluate windfield based on values found in domain
1397        """
1398
1399        from math import pi, cos, sin, sqrt
1400        from Numeric import Float, ones, ArrayType
1401
1402        xmom_update = domain.quantities['xmomentum'].explicit_update
1403        ymom_update = domain.quantities['ymomentum'].explicit_update
1404
[3928]1405        N = len(domain) # number_of_triangles
[3804]1406        t = domain.time
1407
1408        if callable(self.speed):
1409            xc = domain.get_centroid_coordinates()
1410            s_vec = self.speed(t, xc[:,0], xc[:,1])
1411        else:
[4769]1412            # Assume s is a scalar
[3804]1413
1414            try:
1415                s_vec = self.speed * ones(N, Float)
1416            except:
1417                msg = 'Speed must be either callable or a scalar: %s' %self.s
1418                raise msg
1419
1420
1421        if callable(self.phi):
1422            xc = domain.get_centroid_coordinates()
1423            phi_vec = self.phi(t, xc[:,0], xc[:,1])
1424        else:
[4769]1425            # Assume phi is a scalar
[3804]1426
1427            try:
1428                phi_vec = self.phi * ones(N, Float)
1429            except:
1430                msg = 'Angle must be either callable or a scalar: %s' %self.phi
1431                raise msg
1432
1433        assign_windfield_values(xmom_update, ymom_update,
1434                                s_vec, phi_vec, self.const)
1435
1436
1437def assign_windfield_values(xmom_update, ymom_update,
1438                            s_vec, phi_vec, const):
1439    """Python version of assigning wind field to update vectors.
1440    A c version also exists (for speed)
1441    """
1442    from math import pi, cos, sin, sqrt
1443
1444    N = len(s_vec)
1445    for k in range(N):
1446        s = s_vec[k]
1447        phi = phi_vec[k]
1448
[4769]1449        # Convert to radians
[3804]1450        phi = phi*pi/180
1451
[4769]1452        # Compute velocity vector (u, v)
[3804]1453        u = s*cos(phi)
1454        v = s*sin(phi)
1455
[4769]1456        # Compute wind stress
[3804]1457        S = const * sqrt(u**2 + v**2)
1458        xmom_update[k] += S*u
1459        ymom_update[k] += S*v
1460
1461
1462
[5294]1463
1464
1465class General_forcing:
1466    """Class General_forcing - general explicit forcing term for update of quantity
1467   
1468    This is used by Inflow and Rainfall for instance
1469   
1470
1471    General_forcing(quantity_name, rate, center, radius, polygon)
1472
1473    domain:     ANUGA computational domain
1474    quantity_name: Name of quantity to update. It must be a known conserved quantity.
1475    rate [?/s]: Total rate of change over the specified area. 
1476                This parameter can be either a constant or a
1477                function of time. Positive values indicate increases,
1478                negative values indicate decreases.
1479                Rate can be None at initialisation but must be specified
[5436]1480                before forcing term is applied (i.e. simulation has started).
[5294]1481
1482    center [m]: Coordinates at center of flow point
1483    radius [m]: Size of circular area
1484    polygon:    Arbitrary polygon.
1485
1486
1487    Either center, radius or polygon can be specified but not both.
1488    If neither are specified the entire domain gets updated.
1489   
1490    See Inflow or Rainfall for examples of use
1491    """
1492
1493
1494    # FIXME (AnyOne) : Add various methods to allow spatial variations
1495
1496    def __init__(self,
1497                 domain,
1498                 quantity_name,
[5303]1499                 rate=0.0,
[5294]1500                 center=None, radius=None,
1501                 polygon=None,
1502                 verbose=False):
1503                     
[5450]1504        if center is None:
1505            msg = 'I got radius but no center.'       
1506            assert radius is None, msg
1507           
1508        if radius is None:
1509            msg += 'I got center but no radius.'       
1510            assert center is None, msg
1511           
1512           
1513                     
[5570]1514        from math import pi, cos, sin
[5294]1515
1516        self.domain = domain
1517        self.quantity_name = quantity_name
1518        self.rate = rate
1519        self.center = ensure_numeric(center)
1520        self.radius = radius
1521        self.polygon = polygon       
1522        self.verbose = verbose
[5303]1523        self.value = 0.0 # Can be used to remember value at
1524                         # previous timestep in order to obtain rate
[5294]1525
[5570]1526                         
1527        bounding_polygon = domain.get_boundary_polygon()
1528
1529
[5294]1530        # Update area if applicable
[5436]1531        self.exchange_area = None       
[5294]1532        if center is not None and radius is not None:
1533            assert len(center) == 2
1534            msg = 'Polygon cannot be specified when center and radius are'
1535            assert polygon is None, msg
1536
[5436]1537            self.exchange_area = radius**2*pi
[5570]1538
1539            # Check that circle center lies within the mesh.
1540            msg = 'Center %s specified for forcing term did not' %(str(center))
1541            msg += 'fall within the domain boundary.'
1542            assert is_inside_polygon(center, bounding_polygon), msg
1543
1544            # Check that circle periphery lies within the mesh.
1545            N = 100
1546            periphery_points = []
1547            for i in range(N):
1548
1549                theta = 2*pi*i/100
1550               
1551                x = center[0] + radius*cos(theta)
1552                y = center[1] + radius*sin(theta)
1553
1554                periphery_points.append([x,y])
1555
1556
1557            for point in periphery_points:
1558                msg = 'Point %s on periphery for forcing term did not' %(str(point))
1559                msg += ' fall within the domain boundary.'
1560                assert is_inside_polygon(point, bounding_polygon), msg
1561
[5294]1562       
1563        if polygon is not None:
[5436]1564            self.exchange_area = polygon_area(self.polygon)
[5294]1565
[5570]1566            # Check that polygon lies within the mesh.
1567            for point in self.polygon:
1568                msg = 'Point %s in polygon for forcing term did not' %(str(point))
1569                msg += 'fall within the domain boundary.'
1570                assert is_inside_polygon(point, bounding_polygon), msg
1571       
[5294]1572
[5570]1573
1574           
1575
1576
[5294]1577        # Pointer to update vector
1578        self.update = domain.quantities[self.quantity_name].explicit_update           
1579
1580        # Determine indices in flow area
1581        N = len(domain)   
1582        points = domain.get_centroid_coordinates(absolute=True)
1583
[5436]1584        # Calculate indices in exchange area for this forcing term
1585        self.exchange_indices = None
[5294]1586        if self.center is not None and self.radius is not None:
1587            # Inlet is circular
1588           
[5450]1589            inlet_region = 'center=%s, radius=%s' %(self.center, self.radius) 
1590           
[5436]1591            self.exchange_indices = []
[5294]1592            for k in range(N):
1593                x, y = points[k,:] # Centroid
1594                if ((x-self.center[0])**2+(y-self.center[1])**2) < self.radius**2:
[5436]1595                    self.exchange_indices.append(k)
[5294]1596                   
1597        if self.polygon is not None:                   
1598            # Inlet is polygon
[5450]1599           
1600            inlet_region = 'polygon=%s' %(self.polygon) 
1601                       
[5436]1602            self.exchange_indices = inside_polygon(points, self.polygon)
[5294]1603           
1604           
[5450]1605        if self.exchange_indices is not None:
1606            #print inlet_region
1607       
1608            if len(self.exchange_indices) == 0:
1609                msg = 'No triangles have been identified in specified region: %s' %inlet_region
1610                raise Exception, msg
[5294]1611
1612
1613    def __call__(self, domain):
1614        """Apply inflow function at time specified in domain and update stage
1615        """
1616
1617        # Call virtual method allowing local modifications
1618        rate = self.update_rate(domain.get_time())
1619        if rate is None:
1620            msg = 'Attribute rate must be specified in General_forcing'
1621            msg += ' or its descendants before attempting to call it'
1622            raise Exception, msg
1623       
1624
1625        # Now rate is a number
1626        if self.verbose is True:
1627            print 'Rate of %s at time = %.2f = %f' %(self.quantity_name,
1628                                                     domain.get_time(),
1629                                                     rate)
1630
1631
[5436]1632        if self.exchange_indices is None:
[5294]1633            self.update[:] += rate
1634        else:
1635            # Brute force assignment of restricted rate
[5436]1636            for k in self.exchange_indices:
[5294]1637                self.update[k] += rate
1638
1639
1640    def update_rate(self, t):
1641        """Virtual method allowing local modifications by writing an
1642        overriding version in descendant
1643       
1644        """
1645        if callable(self.rate):
1646            rate = self.rate(t)
1647        else:
1648            rate = self.rate
1649
1650        return rate
1651
1652
1653    def get_quantity_values(self):
1654        """Return values for specified quantity restricted to opening
1655        """
[5436]1656        return self.domain.quantities[self.quantity_name].get_values(indices=self.exchange_indices)
[5294]1657   
1658
1659    def set_quantity_values(self, val):
1660        """Set values for specified quantity restricted to opening
1661        """
[5436]1662        self.domain.quantities[self.quantity_name].set_values(val, indices=self.exchange_indices)   
[5294]1663
1664
1665
1666class Rainfall(General_forcing):
[4530]1667    """Class Rainfall - general 'rain over entire domain' forcing term.
1668   
1669    Used for implementing Rainfall over the entire domain.
1670       
1671        Current Limited to only One Gauge..
1672       
1673        Need to add Spatial Varying Capability
1674        (This module came from copying and amending the Inflow Code)
1675   
1676    Rainfall(rain)
[5294]1677
1678    domain   
[4530]1679    rain [mm/s]:  Total rain rate over the specified domain. 
1680                  NOTE: Raingauge Data needs to reflect the time step.
1681                  IE: if Gauge is mm read at a time step, then the input
1682                  here is as mm/(timeStep) so 10mm in 5minutes becomes
1683                  10/(5x60) = 0.0333mm/s.
1684       
1685       
1686                  This parameter can be either a constant or a
1687                  function of time. Positive values indicate inflow,
1688                  negative values indicate outflow.
1689                  (and be used for Infiltration - Write Seperate Module)
1690                  The specified flow will be divided by the area of
1691                  the inflow region and then applied to update the
[5294]1692                  stage quantity.
[5178]1693
1694    polygon: Specifies a polygon to restrict the rainfall.
[4530]1695   
1696    Examples
1697    How to put them in a run File...
1698       
1699    #--------------------------------------------------------------------------
1700    # Setup specialised forcing terms
1701    #--------------------------------------------------------------------------
1702    # This is the new element implemented by Ole and Rudy to allow direct
1703    # input of Inflow in mm/s
[4438]1704
[4530]1705    catchmentrainfall = Rainfall(rain=file_function('Q100_2hr_Rain.tms')) 
1706                        # Note need path to File in String.
1707                        # Else assumed in same directory
1708
1709    domain.forcing_terms.append(catchmentrainfall)
1710    """
1711
[5178]1712   
[4530]1713    def __init__(self,
[5294]1714                 domain,
1715                 rate=0.0,
1716                 center=None, radius=None,
1717                 polygon=None,
1718                 verbose=False):
[4530]1719
[5294]1720        # Converting mm/s to m/s to apply in ANUGA)
1721        if callable(rate):
1722            rain = lambda t: rate(t)/1000.0
1723        else:
1724            rain = rate/1000.0           
[5178]1725           
[5294]1726        General_forcing.__init__(self,
1727                                 domain,
1728                                 'stage',
1729                                 rate=rain,
1730                                 center=center, radius=radius,
1731                                 polygon=polygon,
1732                                 verbose=verbose)
[5178]1733
1734       
[4530]1735
[5178]1736
1737
1738
[5294]1739class Inflow(General_forcing):
[4438]1740    """Class Inflow - general 'rain and drain' forcing term.
1741   
1742    Useful for implementing flows in and out of the domain.
1743   
[5294]1744    Inflow(flow, center, radius, polygon)
1745
1746    domain
[5506]1747    rate [m^3/s]: Total flow rate over the specified area. 
[4438]1748                  This parameter can be either a constant or a
1749                  function of time. Positive values indicate inflow,
1750                  negative values indicate outflow.
1751                  The specified flow will be divided by the area of
[5506]1752                  the inflow region and then applied to update stage.     
[5294]1753    center [m]: Coordinates at center of flow point
1754    radius [m]: Size of circular area
1755    polygon:    Arbitrary polygon.
1756
1757    Either center, radius or polygon must be specified
[4438]1758   
1759    Examples
1760
1761    # Constant drain at 0.003 m^3/s.
1762    # The outflow area is 0.07**2*pi=0.0154 m^2
1763    # This corresponds to a rate of change of 0.003/0.0154 = 0.2 m/s
1764    #                                     
1765    Inflow((0.7, 0.4), 0.07, -0.003)
1766
1767
1768    # Tap turning up to a maximum inflow of 0.0142 m^3/s.
1769    # The inflow area is 0.03**2*pi = 0.00283 m^2
1770    # This corresponds to a rate of change of 0.0142/0.00283 = 5 m/s     
1771    # over the specified area
1772    Inflow((0.5, 0.5), 0.03, lambda t: min(0.01*t, 0.0142))
[4530]1773
[5294]1774
[4530]1775    #--------------------------------------------------------------------------
1776    # Setup specialised forcing terms
1777    #--------------------------------------------------------------------------
1778    # This is the new element implemented by Ole to allow direct input
1779    # of Inflow in m^3/s
1780
1781    hydrograph = Inflow(center=(320, 300), radius=10,
[5506]1782                        rate=file_function('Q/QPMF_Rot_Sub13.tms'))
[4530]1783
1784    domain.forcing_terms.append(hydrograph)
1785   
[4438]1786    """
1787
1788
1789    def __init__(self,
[5294]1790                 domain,
1791                 rate=0.0,
[4438]1792                 center=None, radius=None,
[5294]1793                 polygon=None,
1794                 verbose=False):                 
[4438]1795
1796
[5294]1797        #msg = 'Class Inflow must have either center & radius or a polygon specified.'
1798        #assert center is not None and radius is not None or\
1799        #       polygon is not None, msg
1800
1801
1802        # Create object first to make area is available
1803        General_forcing.__init__(self,
1804                                 domain,
1805                                 'stage',
1806                                 rate=rate,
1807                                 center=center, radius=radius,
1808                                 polygon=polygon,
1809                                 verbose=verbose)
1810
1811    def update_rate(self, t):
1812        """Virtual method allowing local modifications by writing an
1813        overriding version in descendant
1814
1815        This one converts m^3/s to m/s which can be added directly to 'stage' in ANUGA
1816        """
1817
[4438]1818       
[5294]1819       
1820        if callable(self.rate):
[5436]1821            _rate = self.rate(t)/self.exchange_area
[4438]1822        else:
[5436]1823            _rate = self.rate/self.exchange_area
[4438]1824
[5294]1825        return _rate
[4438]1826
1827
1828
1829
[4733]1830#------------------
1831# Initialise module
1832#------------------
[3804]1833
[4769]1834
[3804]1835from anuga.utilities import compile
1836if compile.can_use_C_extension('shallow_water_ext.c'):
[4769]1837    # Underlying C implementations can be accessed
[3804]1838
1839    from shallow_water_ext import rotate, assign_windfield_values
[4769]1840else:
1841    msg = 'C implementations could not be accessed by %s.\n ' %__file__
1842    msg += 'Make sure compile_all.py has been run as described in '
1843    msg += 'the ANUGA installation guide.'
1844    raise Exception, msg
[3804]1845
1846
[4733]1847# Optimisation with psyco
[3804]1848from anuga.config import use_psyco
1849if use_psyco:
1850    try:
1851        import psyco
1852    except:
1853        import os
1854        if os.name == 'posix' and os.uname()[4] == 'x86_64':
1855            pass
1856            #Psyco isn't supported on 64 bit systems, but it doesn't matter
1857        else:
1858            msg = 'WARNING: psyco (speedup) could not import'+\
1859                  ', you may want to consider installing it'
1860            print msg
1861    else:
1862        psyco.bind(Domain.distribute_to_vertices_and_edges)
1863        psyco.bind(Domain.compute_fluxes)
1864
1865if __name__ == "__main__":
1866    pass
1867
[4733]1868
Note: See TracBrowser for help on using the repository browser.