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

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

Rolled back changes from changeset:5297 regarding ticket:262.
All tests and validations passes.

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