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

Last change on this file since 5306 was 5306, checked in by steve, 16 years ago

Working version of wave.py for Ole to work on.

  • 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: 5306 $)
70ModifiedBy:
71    $Author: steve $
72    $Date: 2008-05-10 22:40:01 +0000 (Sat, 10 May 2008) $
73
74"""
75
76# Subversion keywords:
77#
78# $LastChangedDate: 2008-05-10 22:40:01 +0000 (Sat, 10 May 2008) $
79# $LastChangedRevision: 5306 $
80# $LastChangedBy: steve $
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              int(domain.use_centroid_velocities))
736
737
738def distribute_using_vertex_limiter(domain):
739    """Distribution from centroids to vertices specific to the
740    shallow water wave
741    equation.
742
743    It will ensure that h (w-z) is always non-negative even in the
744    presence of steep bed-slopes by taking a weighted average between shallow
745    and deep cases.
746
747    In addition, all conserved quantities get distributed as per either a
748    constant (order==1) or a piecewise linear function (order==2).
749
750    FIXME: more explanation about removal of artificial variability etc
751
752    Precondition:
753      All quantities defined at centroids and bed elevation defined at
754      vertices.
755
756    Postcondition
757      Conserved quantities defined at vertices
758
759    """
760
761   
762
763    # Remove very thin layers of water
764    protect_against_infinitesimal_and_negative_heights(domain)
765
766    # Extrapolate all conserved quantities
767    if domain.optimised_gradient_limiter:
768        # MH090605 if second order,
769        # perform the extrapolation and limiting on
770        # all of the conserved quantities
771
772        if (domain._order_ == 1):
773            for name in domain.conserved_quantities:
774                Q = domain.quantities[name]
775                Q.extrapolate_first_order()
776        elif domain._order_ == 2:
777            domain.extrapolate_second_order_sw()
778        else:
779            raise 'Unknown order'
780    else:
781        # Old code:
782        for name in domain.conserved_quantities:
783            Q = domain.quantities[name]
784
785            if domain._order_ == 1:
786                Q.extrapolate_first_order()
787            elif domain._order_ == 2:
788                Q.extrapolate_second_order_and_limit_by_vertex()
789            else:
790                raise 'Unknown order'
791
792
793    # Take bed elevation into account when water heights are small
794    balance_deep_and_shallow(domain)
795
796    # Compute edge values by interpolation
797    for name in domain.conserved_quantities:
798        Q = domain.quantities[name]
799        Q.interpolate_from_vertices_to_edges()
800
801
802
803def distribute_using_edge_limiter(domain):
804    """Distribution from centroids to edges specific to the
805    shallow water wave
806    equation.
807
808    It will ensure that h (w-z) is always non-negative even in the
809    presence of steep bed-slopes by taking a weighted average between shallow
810    and deep cases.
811
812    In addition, all conserved quantities get distributed as per either a
813    constant (order==1) or a piecewise linear function (order==2).
814
815
816    Precondition:
817      All quantities defined at centroids and bed elevation defined at
818      vertices.
819
820    Postcondition
821      Conserved quantities defined at vertices
822
823    """
824
825    # Remove very thin layers of water
826    protect_against_infinitesimal_and_negative_heights(domain)
827
828
829    for name in domain.conserved_quantities:
830        Q = domain.quantities[name]
831        if domain._order_ == 1:
832            Q.extrapolate_first_order()
833        elif domain._order_ == 2:
834            Q.extrapolate_second_order_and_limit_by_edge()
835        else:
836            raise 'Unknown order'
837
838    balance_deep_and_shallow(domain)
839
840    # Compute edge values by interpolation
841    for name in domain.conserved_quantities:
842        Q = domain.quantities[name]
843        Q.interpolate_from_vertices_to_edges()
844
845
846def protect_against_infinitesimal_and_negative_heights(domain):
847    """Protect against infinitesimal heights and associated high velocities
848    """
849
850    # Shortcuts
851    wc = domain.quantities['stage'].centroid_values
852    zc = domain.quantities['elevation'].centroid_values
853    xmomc = domain.quantities['xmomentum'].centroid_values
854    ymomc = domain.quantities['ymomentum'].centroid_values
855
856    from shallow_water_ext import protect
857
858    protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
859            domain.epsilon, wc, zc, xmomc, ymomc)
860
861
862def h_limiter(domain):
863    """Limit slopes for each volume to eliminate artificial variance
864    introduced by e.g. second order extrapolator
865
866    limit on h = w-z
867
868    This limiter depends on two quantities (w,z) so it resides within
869    this module rather than within quantity.py
870
871    Wrapper for c-extension
872    """
873
874    N = len(domain) # number_of_triangles
875    beta_h = domain.beta_h
876
877    # Shortcuts
878    wc = domain.quantities['stage'].centroid_values
879    zc = domain.quantities['elevation'].centroid_values
880    hc = wc - zc
881
882    wv = domain.quantities['stage'].vertex_values
883    zv = domain.quantities['elevation'].vertex_values
884    hv = wv - zv
885
886    #Call C-extension
887    from shallow_water_ext import h_limiter_sw
888    hvbar = h_limiter_sw(domain, hc, hv)
889
890    return hvbar
891
892
893def balance_deep_and_shallow(domain):
894    """Compute linear combination between stage as computed by
895    gradient-limiters limiting using w, and stage computed by
896    gradient-limiters limiting using h (h-limiter).
897    The former takes precedence when heights are large compared to the
898    bed slope while the latter takes precedence when heights are
899    relatively small.  Anything in between is computed as a balanced
900    linear combination in order to avoid numerical disturbances which
901    would otherwise appear as a result of hard switching between
902    modes.
903
904    Wrapper for C implementation
905    """
906
907    # FIXME (Ole): I reckon this can be simplified significantly:
908    #
909    # Always use beta_h == 0, and phase it out.
910    # Compute hc and hv in the c-code
911    # Omit updating xmomv
912    #
913    from shallow_water_ext import balance_deep_and_shallow as balance_deep_and_shallow_c
914
915
916    #print 'calling balance depth and shallow'
917    # Shortcuts
918    wc = domain.quantities['stage'].centroid_values
919    zc = domain.quantities['elevation'].centroid_values
920
921    wv = domain.quantities['stage'].vertex_values
922    zv = domain.quantities['elevation'].vertex_values
923
924    # Momentums at centroids
925    xmomc = domain.quantities['xmomentum'].centroid_values
926    ymomc = domain.quantities['ymomentum'].centroid_values
927
928    # Momentums at vertices
929    xmomv = domain.quantities['xmomentum'].vertex_values
930    ymomv = domain.quantities['ymomentum'].vertex_values
931
932    # Limit h
933    if domain.beta_h > 0:
934        hvbar = h_limiter(domain)
935       
936        balance_deep_and_shallow_c(domain, domain.beta_h,
937                                   wc, zc, wv, zv, hvbar,
938                                   xmomc, ymomc, xmomv, ymomv)       
939    else:
940        # print 'Using first order h-limiter'
941        # FIXME: Pass wc in for now - it will be ignored.
942       
943        # This is how one would make a first order h_limited value
944        # as in the old balancer (pre 17 Feb 2005):
945        #  If we wish to hard wire this, one should modify the C-code
946        # from Numeric import zeros, Float
947        # hvbar = zeros( (len(wc), 3), Float)
948        # for i in range(3):
949        #     hvbar[:,i] = wc[:] - zc[:]
950
951        balance_deep_and_shallow_c(domain, domain.beta_h,
952                                   wc, zc, wv, zv, wc, 
953                                   xmomc, ymomc, xmomv, ymomv)
954
955
956
957
958#------------------------------------------------------------------
959# Boundary conditions - specific to the shallow water wave equation
960#------------------------------------------------------------------
961class Reflective_boundary(Boundary):
962    """Reflective boundary returns same conserved quantities as
963    those present in its neighbour volume but reflected.
964
965    This class is specific to the shallow water equation as it
966    works with the momentum quantities assumed to be the second
967    and third conserved quantities.
968    """
969
970    def __init__(self, domain = None):
971        Boundary.__init__(self)
972
973        if domain is None:
974            msg = 'Domain must be specified for reflective boundary'
975            raise msg
976
977        # Handy shorthands
978        self.stage   = domain.quantities['stage'].edge_values
979        self.xmom    = domain.quantities['xmomentum'].edge_values
980        self.ymom    = domain.quantities['ymomentum'].edge_values
981        self.normals = domain.normals
982
983        self.conserved_quantities = zeros(3, Float)
984
985    def __repr__(self):
986        return 'Reflective_boundary'
987
988
989    def evaluate(self, vol_id, edge_id):
990        """Reflective boundaries reverses the outward momentum
991        of the volume they serve.
992        """
993
994        q = self.conserved_quantities
995        q[0] = self.stage[vol_id, edge_id]
996        q[1] = self.xmom[vol_id, edge_id]
997        q[2] = self.ymom[vol_id, edge_id]
998
999        normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
1000
1001
1002        r = rotate(q, normal, direction = 1)
1003        r[1] = -r[1]
1004        q = rotate(r, normal, direction = -1)
1005
1006        return q
1007
1008
1009
1010class Transmissive_Momentum_Set_Stage_boundary(Boundary):
1011    """Returns same momentum conserved quantities as
1012    those present in its neighbour volume.
1013    Sets stage by specifying a function f of time which may either be a
1014    vector function or a scalar function
1015
1016    Example:
1017
1018    def waveform(t):
1019        return sea_level + normalized_amplitude/cosh(t-25)**2
1020
1021    Bts = Transmissive_Momentum_Set_Stage_boundary(domain, waveform)
1022   
1023
1024    Underlying domain must be specified when boundary is instantiated
1025    """
1026
1027    def __init__(self, domain = None, function=None):
1028        Boundary.__init__(self)
1029
1030        if domain is None:
1031            msg = 'Domain must be specified for this type boundary'
1032            raise msg
1033
1034        if function is None:
1035            msg = 'Function must be specified for this type boundary'
1036            raise msg
1037
1038        self.domain   = domain
1039        self.function = function
1040
1041    def __repr__(self):
1042        return 'Transmissive_Momentum_Set_Stage_boundary(%s)' %self.domain
1043
1044    def evaluate(self, vol_id, edge_id):
1045        """Transmissive Momentum Set Stage boundaries return the edge momentum
1046        values of the volume they serve.
1047        """
1048
1049        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
1050
1051
1052        t = self.domain.time
1053
1054        if hasattr(self.function, 'time'):
1055            # Roll boundary over if time exceeds           
1056            while t > self.function.time[-1]:
1057                msg = 'WARNING: domain time %.2f has exceeded' %t
1058                msg += 'time provided in '
1059                msg += 'transmissive_momentum_set_stage boundary object.\n'
1060                msg += 'I will continue, reusing the object from t==0'
1061                print msg
1062                t -= self.function.time[-1]
1063
1064
1065        value = self.function(t)
1066        try:
1067            x = float(value)
1068        except:
1069            x = float(value[0])
1070           
1071        q[0] = x
1072        return q
1073
1074
1075        # FIXME: Consider this (taken from File_boundary) to allow
1076        # spatial variation
1077        # if vol_id is not None and edge_id is not None:
1078        #     i = self.boundary_indices[ vol_id, edge_id ]
1079        #     return self.F(t, point_id = i)
1080        # else:
1081        #     return self.F(t)
1082
1083
1084
1085class Dirichlet_Discharge_boundary(Boundary):
1086    """
1087    Sets stage (stage0)
1088    Sets momentum (wh0) in the inward normal direction.
1089
1090    Underlying domain must be specified when boundary is instantiated
1091    """
1092
1093    def __init__(self, domain = None, stage0=None, wh0=None):
1094        Boundary.__init__(self)
1095
1096        if domain is None:
1097            msg = 'Domain must be specified for this type boundary'
1098            raise msg
1099
1100        if stage0 is None:
1101            raise 'set stage'
1102
1103        if wh0 is None:
1104            wh0 = 0.0
1105
1106        self.domain   = domain
1107        self.stage0  = stage0
1108        self.wh0 = wh0
1109
1110    def __repr__(self):
1111        return 'Dirichlet_Discharge_boundary(%s)' %self.domain
1112
1113    def evaluate(self, vol_id, edge_id):
1114        """Set discharge in the (inward) normal direction
1115        """
1116
1117        normal = self.domain.get_normal(vol_id,edge_id)
1118        q = [self.stage0, -self.wh0*normal[0], -self.wh0*normal[1]]
1119        return q
1120
1121
1122        # FIXME: Consider this (taken from File_boundary) to allow
1123        # spatial variation
1124        # if vol_id is not None and edge_id is not None:
1125        #     i = self.boundary_indices[ vol_id, edge_id ]
1126        #     return self.F(t, point_id = i)
1127        # else:
1128        #     return self.F(t)
1129
1130
1131class Field_boundary(Boundary):
1132    """Set boundary from given field represented in an sww file containing values
1133    for stage, xmomentum and ymomentum.
1134    Optionally, the user can specify mean_stage to offset the stage provided in the
1135    sww file.
1136
1137    This function is a thin wrapper around the generic File_boundary. The
1138    difference between the file_boundary and field_boundary is only that the
1139    field_boundary will allow you to change the level of the stage height when
1140    you read in the boundary condition. This is very useful when running
1141    different tide heights in the same area as you need only to convert one
1142    boundary condition to a SWW file, ideally for tide height of 0 m
1143    (saving disk space). Then you can use field_boundary to read this SWW file
1144    and change the stage height (tide) on the fly depending on the scenario.
1145   
1146    """
1147
1148
1149    def __init__(self, filename, domain,
1150                 mean_stage=0.0,
1151                 time_thinning=1, 
1152                 use_cache=False,
1153                 verbose=False):
1154        """Constructor
1155
1156        filename: Name of sww file
1157        domain: pointer to shallow water domain for which the boundary applies
1158        mean_stage: The mean water level which will be added to stage derived
1159                    from the sww file
1160        time_thinning: Will set how many time steps from the sww file read in
1161                       will be interpolated to the boundary. For example if
1162                       the sww file has 1 second time steps and is 24 hours
1163                       in length it has 86400 time steps. If you set
1164                       time_thinning to 1 it will read all these steps.
1165                       If you set it to 100 it will read every 100th step eg
1166                       only 864 step. This parameter is very useful to increase
1167                       the speed of a model run that you are setting up
1168                       and testing.
1169        use_cache:
1170        verbose:
1171       
1172        """
1173
1174        # Create generic file_boundary object
1175        self.file_boundary = File_boundary(filename, domain,
1176                                           time_thinning=time_thinning,
1177                                           use_cache=use_cache,
1178                                           verbose=verbose)
1179       
1180        # Record information from File_boundary
1181        self.F = self.file_boundary.F
1182        self.domain = self.file_boundary.domain
1183       
1184        # Record mean stage
1185        self.mean_stage = mean_stage
1186
1187
1188    def __repr__(self):
1189        return 'Field boundary'
1190
1191
1192    def evaluate(self, vol_id=None, edge_id=None):
1193        """Return linearly interpolated values based on domain.time
1194
1195        vol_id and edge_id are ignored
1196        """
1197
1198        # Evaluate file boundary
1199        q = self.file_boundary.evaluate(vol_id, edge_id)
1200
1201        # Adjust stage
1202        for j, name in enumerate(self.domain.conserved_quantities):
1203            if name == 'stage':
1204                q[j] += self.mean_stage
1205        return q
1206
1207   
1208
1209#-----------------------
1210# Standard forcing terms
1211#-----------------------
1212
1213def gravity(domain):
1214    """Apply gravitational pull in the presence of bed slope
1215    Wrapper calls underlying C implementation
1216    """
1217
1218    xmom = domain.quantities['xmomentum'].explicit_update
1219    ymom = domain.quantities['ymomentum'].explicit_update
1220
1221    stage = domain.quantities['stage']
1222    elevation = domain.quantities['elevation']
1223
1224    h = stage.centroid_values - elevation.centroid_values
1225    z = elevation.vertex_values
1226
1227    x = domain.get_vertex_coordinates()
1228    g = domain.g
1229   
1230
1231    from shallow_water_ext import gravity as gravity_c
1232    gravity_c(g, h, z, x, xmom, ymom) #, 1.0e-6)
1233
1234
1235
1236def manning_friction_implicit(domain):
1237    """Apply (Manning) friction to water momentum   
1238    Wrapper for c version
1239    """
1240
1241
1242    #print 'Implicit friction'
1243
1244    xmom = domain.quantities['xmomentum']
1245    ymom = domain.quantities['ymomentum']
1246
1247    w = domain.quantities['stage'].centroid_values
1248    z = domain.quantities['elevation'].centroid_values
1249
1250    uh = xmom.centroid_values
1251    vh = ymom.centroid_values
1252    eta = domain.quantities['friction'].centroid_values
1253
1254    xmom_update = xmom.semi_implicit_update
1255    ymom_update = ymom.semi_implicit_update
1256
1257    N = len(domain)
1258    eps = domain.minimum_allowed_height
1259    g = domain.g
1260
1261    from shallow_water_ext import manning_friction as manning_friction_c
1262    manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
1263
1264
1265def manning_friction_explicit(domain):
1266    """Apply (Manning) friction to water momentum   
1267    Wrapper for c version
1268    """
1269
1270    # print 'Explicit friction'
1271
1272    xmom = domain.quantities['xmomentum']
1273    ymom = domain.quantities['ymomentum']
1274
1275    w = domain.quantities['stage'].centroid_values
1276    z = domain.quantities['elevation'].centroid_values
1277
1278    uh = xmom.centroid_values
1279    vh = ymom.centroid_values
1280    eta = domain.quantities['friction'].centroid_values
1281
1282    xmom_update = xmom.explicit_update
1283    ymom_update = ymom.explicit_update
1284
1285    N = len(domain)
1286    eps = domain.minimum_allowed_height
1287    g = domain.g
1288
1289    from shallow_water_ext import manning_friction as manning_friction_c
1290    manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
1291
1292
1293# FIXME (Ole): This was implemented for use with one of the analytical solutions (Sampson?)
1294# Is it still needed (30 Oct 2007)?
1295def linear_friction(domain):
1296    """Apply linear friction to water momentum
1297
1298    Assumes quantity: 'linear_friction' to be present
1299    """
1300
1301    from math import sqrt
1302
1303    w = domain.quantities['stage'].centroid_values
1304    z = domain.quantities['elevation'].centroid_values
1305    h = w-z
1306
1307    uh = domain.quantities['xmomentum'].centroid_values
1308    vh = domain.quantities['ymomentum'].centroid_values
1309    tau = domain.quantities['linear_friction'].centroid_values
1310
1311    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1312    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1313
1314    N = len(domain) # number_of_triangles
1315    eps = domain.minimum_allowed_height
1316    g = domain.g #Not necessary? Why was this added?
1317
1318    for k in range(N):
1319        if tau[k] >= eps:
1320            if h[k] >= eps:
1321                S = -tau[k]/h[k]
1322
1323                #Update momentum
1324                xmom_update[k] += S*uh[k]
1325                ymom_update[k] += S*vh[k]
1326
1327
1328
1329#---------------------------------
1330# Experimental auxiliary functions
1331#---------------------------------
1332def check_forcefield(f):
1333    """Check that f is either
1334    1: a callable object f(t,x,y), where x and y are vectors
1335       and that it returns an array or a list of same length
1336       as x and y
1337    2: a scalar
1338    """
1339
1340    if callable(f):
1341        N = 3
1342        x = ones(3, Float)
1343        y = ones(3, Float)
1344        try:
1345            q = f(1.0, x=x, y=y)
1346        except Exception, e:
1347            msg = 'Function %s could not be executed:\n%s' %(f, e)
1348            # FIXME: Reconsider this semantics
1349            raise msg
1350
1351        try:
1352            q = array(q).astype(Float)
1353        except:
1354            msg = 'Return value from vector function %s could ' %f
1355            msg += 'not be converted into a Numeric array of floats.\n'
1356            msg += 'Specified function should return either list or array.'
1357            raise msg
1358
1359        # Is this really what we want?
1360        msg = 'Return vector from function %s ' %f
1361        msg += 'must have same lenght as input vectors'
1362        assert len(q) == N, msg
1363
1364    else:
1365        try:
1366            f = float(f)
1367        except:
1368            msg = 'Force field %s must be either a scalar' %f
1369            msg += ' or a vector function'
1370            raise Exception(msg)
1371    return f
1372
1373
1374class Wind_stress:
1375    """Apply wind stress to water momentum in terms of
1376    wind speed [m/s] and wind direction [degrees]
1377    """
1378
1379    def __init__(self, *args, **kwargs):
1380        """Initialise windfield from wind speed s [m/s]
1381        and wind direction phi [degrees]
1382
1383        Inputs v and phi can be either scalars or Python functions, e.g.
1384
1385        W = Wind_stress(10, 178)
1386
1387        #FIXME - 'normal' degrees are assumed for now, i.e. the
1388        vector (1,0) has zero degrees.
1389        We may need to convert from 'compass' degrees later on and also
1390        map from True north to grid north.
1391
1392        Arguments can also be Python functions of t,x,y as in
1393
1394        def speed(t,x,y):
1395            ...
1396            return s
1397
1398        def angle(t,x,y):
1399            ...
1400            return phi
1401
1402        where x and y are vectors.
1403
1404        and then pass the functions in
1405
1406        W = Wind_stress(speed, angle)
1407
1408        The instantiated object W can be appended to the list of
1409        forcing_terms as in
1410
1411        Alternatively, one vector valued function for (speed, angle)
1412        can be applied, providing both quantities simultaneously.
1413        As in
1414        W = Wind_stress(F), where returns (speed, angle) for each t.
1415
1416        domain.forcing_terms.append(W)
1417        """
1418
1419        from anuga.config import rho_a, rho_w, eta_w
1420        from Numeric import array, Float
1421
1422        if len(args) == 2:
1423            s = args[0]
1424            phi = args[1]
1425        elif len(args) == 1:
1426            # Assume vector function returning (s, phi)(t,x,y)
1427            vector_function = args[0]
1428            s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
1429            phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
1430        else:
1431           # Assume info is in 2 keyword arguments
1432
1433           if len(kwargs) == 2:
1434               s = kwargs['s']
1435               phi = kwargs['phi']
1436           else:
1437               raise 'Assumes two keyword arguments: s=..., phi=....'
1438
1439        self.speed = check_forcefield(s)
1440        self.phi = check_forcefield(phi)
1441
1442        self.const = eta_w*rho_a/rho_w
1443
1444
1445    def __call__(self, domain):
1446        """Evaluate windfield based on values found in domain
1447        """
1448
1449        from math import pi, cos, sin, sqrt
1450        from Numeric import Float, ones, ArrayType
1451
1452        xmom_update = domain.quantities['xmomentum'].explicit_update
1453        ymom_update = domain.quantities['ymomentum'].explicit_update
1454
1455        N = len(domain) # number_of_triangles
1456        t = domain.time
1457
1458        if callable(self.speed):
1459            xc = domain.get_centroid_coordinates()
1460            s_vec = self.speed(t, xc[:,0], xc[:,1])
1461        else:
1462            # Assume s is a scalar
1463
1464            try:
1465                s_vec = self.speed * ones(N, Float)
1466            except:
1467                msg = 'Speed must be either callable or a scalar: %s' %self.s
1468                raise msg
1469
1470
1471        if callable(self.phi):
1472            xc = domain.get_centroid_coordinates()
1473            phi_vec = self.phi(t, xc[:,0], xc[:,1])
1474        else:
1475            # Assume phi is a scalar
1476
1477            try:
1478                phi_vec = self.phi * ones(N, Float)
1479            except:
1480                msg = 'Angle must be either callable or a scalar: %s' %self.phi
1481                raise msg
1482
1483        assign_windfield_values(xmom_update, ymom_update,
1484                                s_vec, phi_vec, self.const)
1485
1486
1487def assign_windfield_values(xmom_update, ymom_update,
1488                            s_vec, phi_vec, const):
1489    """Python version of assigning wind field to update vectors.
1490    A c version also exists (for speed)
1491    """
1492    from math import pi, cos, sin, sqrt
1493
1494    N = len(s_vec)
1495    for k in range(N):
1496        s = s_vec[k]
1497        phi = phi_vec[k]
1498
1499        # Convert to radians
1500        phi = phi*pi/180
1501
1502        # Compute velocity vector (u, v)
1503        u = s*cos(phi)
1504        v = s*sin(phi)
1505
1506        # Compute wind stress
1507        S = const * sqrt(u**2 + v**2)
1508        xmom_update[k] += S*u
1509        ymom_update[k] += S*v
1510
1511
1512
1513
1514
1515class General_forcing:
1516    """Class General_forcing - general explicit forcing term for update of quantity
1517   
1518    This is used by Inflow and Rainfall for instance
1519   
1520
1521    General_forcing(quantity_name, rate, center, radius, polygon)
1522
1523    domain:     ANUGA computational domain
1524    quantity_name: Name of quantity to update. It must be a known conserved quantity.
1525    rate [?/s]: Total rate of change over the specified area. 
1526                This parameter can be either a constant or a
1527                function of time. Positive values indicate increases,
1528                negative values indicate decreases.
1529                Rate can be None at initialisation but must be specified
1530                before forcting term is applied (i.e. simulation has started).
1531
1532    center [m]: Coordinates at center of flow point
1533    radius [m]: Size of circular area
1534    polygon:    Arbitrary polygon.
1535
1536
1537    Either center, radius or polygon can be specified but not both.
1538    If neither are specified the entire domain gets updated.
1539   
1540    See Inflow or Rainfall for examples of use
1541    """
1542
1543
1544    # FIXME (AnyOne) : Add various methods to allow spatial variations
1545
1546    def __init__(self,
1547                 domain,
1548                 quantity_name,
1549                 rate=0.0,
1550                 center=None, radius=None,
1551                 polygon=None,
1552                 verbose=False):
1553                     
1554
1555        from math import pi
1556
1557        self.domain = domain
1558        self.quantity_name = quantity_name
1559        self.rate = rate
1560        self.center = ensure_numeric(center)
1561        self.radius = radius
1562        self.polygon = polygon       
1563        self.verbose = verbose
1564        self.value = 0.0 # Can be used to remember value at
1565                         # previous timestep in order to obtain rate
1566
1567        # Update area if applicable
1568        self.area = None       
1569        if center is not None and radius is not None:
1570            assert len(center) == 2
1571            msg = 'Polygon cannot be specified when center and radius are'
1572            assert polygon is None, msg
1573
1574            self.area = radius**2*pi
1575       
1576        if polygon is not None:
1577            self.area = polygon_area(self.polygon)
1578
1579
1580        # Pointer to update vector
1581        self.update = domain.quantities[self.quantity_name].explicit_update           
1582
1583        # Determine indices in flow area
1584        N = len(domain)   
1585        points = domain.get_centroid_coordinates(absolute=True)
1586
1587        self.indices = None
1588        if self.center is not None and self.radius is not None:
1589            # Inlet is circular
1590           
1591            self.indices = []
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:
1595                    self.indices.append(k)
1596                   
1597        if self.polygon is not None:                   
1598            # Inlet is polygon
1599            self.indices = inside_polygon(points, self.polygon)
1600           
1601
1602           
1603
1604
1605    def __call__(self, domain):
1606        """Apply inflow function at time specified in domain and update stage
1607        """
1608
1609        # Call virtual method allowing local modifications
1610        rate = self.update_rate(domain.get_time())
1611        if rate is None:
1612            msg = 'Attribute rate must be specified in General_forcing'
1613            msg += ' or its descendants before attempting to call it'
1614            raise Exception, msg
1615       
1616
1617        # Now rate is a number
1618        if self.verbose is True:
1619            print 'Rate of %s at time = %.2f = %f' %(self.quantity_name,
1620                                                     domain.get_time(),
1621                                                     rate)
1622
1623
1624        if self.indices is None:
1625            self.update[:] += rate
1626        else:
1627            # Brute force assignment of restricted rate
1628            for k in self.indices:
1629                self.update[k] += rate
1630
1631
1632    def update_rate(self, t):
1633        """Virtual method allowing local modifications by writing an
1634        overriding version in descendant
1635       
1636        """
1637        if callable(self.rate):
1638            rate = self.rate(t)
1639        else:
1640            rate = self.rate
1641
1642        return rate
1643
1644
1645    def get_quantity_values(self):
1646        """Return values for specified quantity restricted to opening
1647        """
1648        return self.domain.quantities[self.quantity_name].get_values(indices=self.indices)
1649   
1650
1651    def set_quantity_values(self, val):
1652        """Set values for specified quantity restricted to opening
1653        """
1654        self.domain.quantities[self.quantity_name].set_values(val, indices=self.indices)   
1655
1656
1657
1658class Rainfall(General_forcing):
1659    """Class Rainfall - general 'rain over entire domain' forcing term.
1660   
1661    Used for implementing Rainfall over the entire domain.
1662       
1663        Current Limited to only One Gauge..
1664       
1665        Need to add Spatial Varying Capability
1666        (This module came from copying and amending the Inflow Code)
1667   
1668    Rainfall(rain)
1669
1670    domain   
1671    rain [mm/s]:  Total rain rate over the specified domain. 
1672                  NOTE: Raingauge Data needs to reflect the time step.
1673                  IE: if Gauge is mm read at a time step, then the input
1674                  here is as mm/(timeStep) so 10mm in 5minutes becomes
1675                  10/(5x60) = 0.0333mm/s.
1676       
1677       
1678                  This parameter can be either a constant or a
1679                  function of time. Positive values indicate inflow,
1680                  negative values indicate outflow.
1681                  (and be used for Infiltration - Write Seperate Module)
1682                  The specified flow will be divided by the area of
1683                  the inflow region and then applied to update the
1684                  stage quantity.
1685
1686    polygon: Specifies a polygon to restrict the rainfall.
1687   
1688    Examples
1689    How to put them in a run File...
1690       
1691    #--------------------------------------------------------------------------
1692    # Setup specialised forcing terms
1693    #--------------------------------------------------------------------------
1694    # This is the new element implemented by Ole and Rudy to allow direct
1695    # input of Inflow in mm/s
1696
1697    catchmentrainfall = Rainfall(rain=file_function('Q100_2hr_Rain.tms')) 
1698                        # Note need path to File in String.
1699                        # Else assumed in same directory
1700
1701    domain.forcing_terms.append(catchmentrainfall)
1702    """
1703
1704   
1705    def __init__(self,
1706                 domain,
1707                 rate=0.0,
1708                 center=None, radius=None,
1709                 polygon=None,
1710                 verbose=False):
1711
1712        # Converting mm/s to m/s to apply in ANUGA)
1713        if callable(rate):
1714            rain = lambda t: rate(t)/1000.0
1715        else:
1716            rain = rate/1000.0           
1717           
1718        General_forcing.__init__(self,
1719                                 domain,
1720                                 'stage',
1721                                 rate=rain,
1722                                 center=center, radius=radius,
1723                                 polygon=polygon,
1724                                 verbose=verbose)
1725
1726       
1727
1728
1729
1730
1731class Inflow(General_forcing):
1732    """Class Inflow - general 'rain and drain' forcing term.
1733   
1734    Useful for implementing flows in and out of the domain.
1735   
1736    Inflow(flow, center, radius, polygon)
1737
1738    domain
1739    flow [m^3/s]: Total flow rate over the specified area. 
1740                  This parameter can be either a constant or a
1741                  function of time. Positive values indicate inflow,
1742                  negative values indicate outflow.
1743                  The specified flow will be divided by the area of
1744                  the inflow region and then applied to update the
1745                  quantity in question.     
1746    center [m]: Coordinates at center of flow point
1747    radius [m]: Size of circular area
1748    polygon:    Arbitrary polygon.
1749
1750    Either center, radius or polygon must be specified
1751   
1752    Examples
1753
1754    # Constant drain at 0.003 m^3/s.
1755    # The outflow area is 0.07**2*pi=0.0154 m^2
1756    # This corresponds to a rate of change of 0.003/0.0154 = 0.2 m/s
1757    #                                     
1758    Inflow((0.7, 0.4), 0.07, -0.003)
1759
1760
1761    # Tap turning up to a maximum inflow of 0.0142 m^3/s.
1762    # The inflow area is 0.03**2*pi = 0.00283 m^2
1763    # This corresponds to a rate of change of 0.0142/0.00283 = 5 m/s     
1764    # over the specified area
1765    Inflow((0.5, 0.5), 0.03, lambda t: min(0.01*t, 0.0142))
1766
1767
1768    #--------------------------------------------------------------------------
1769    # Setup specialised forcing terms
1770    #--------------------------------------------------------------------------
1771    # This is the new element implemented by Ole to allow direct input
1772    # of Inflow in m^3/s
1773
1774    hydrograph = Inflow(center=(320, 300), radius=10,
1775                        flow=file_function('Q/QPMF_Rot_Sub13.tms'))
1776
1777    domain.forcing_terms.append(hydrograph)
1778   
1779    """
1780
1781
1782    def __init__(self,
1783                 domain,
1784                 rate=0.0,
1785                 center=None, radius=None,
1786                 polygon=None,
1787                 verbose=False):                 
1788
1789
1790        #msg = 'Class Inflow must have either center & radius or a polygon specified.'
1791        #assert center is not None and radius is not None or\
1792        #       polygon is not None, msg
1793
1794
1795        # Create object first to make area is available
1796        General_forcing.__init__(self,
1797                                 domain,
1798                                 'stage',
1799                                 rate=rate,
1800                                 center=center, radius=radius,
1801                                 polygon=polygon,
1802                                 verbose=verbose)
1803
1804    def update_rate(self, t):
1805        """Virtual method allowing local modifications by writing an
1806        overriding version in descendant
1807
1808        This one converts m^3/s to m/s which can be added directly to 'stage' in ANUGA
1809        """
1810
1811       
1812       
1813        if callable(self.rate):
1814            _rate = self.rate(t)/self.area
1815        else:
1816            _rate = self.rate/self.area
1817
1818        return _rate
1819
1820
1821
1822
1823#------------------
1824# Initialise module
1825#------------------
1826
1827
1828from anuga.utilities import compile
1829if compile.can_use_C_extension('shallow_water_ext.c'):
1830    # Underlying C implementations can be accessed
1831
1832    from shallow_water_ext import rotate, assign_windfield_values
1833else:
1834    msg = 'C implementations could not be accessed by %s.\n ' %__file__
1835    msg += 'Make sure compile_all.py has been run as described in '
1836    msg += 'the ANUGA installation guide.'
1837    raise Exception, msg
1838
1839
1840# Optimisation with psyco
1841from anuga.config import use_psyco
1842if use_psyco:
1843    try:
1844        import psyco
1845    except:
1846        import os
1847        if os.name == 'posix' and os.uname()[4] == 'x86_64':
1848            pass
1849            #Psyco isn't supported on 64 bit systems, but it doesn't matter
1850        else:
1851            msg = 'WARNING: psyco (speedup) could not import'+\
1852                  ', you may want to consider installing it'
1853            print msg
1854    else:
1855        psyco.bind(Domain.distribute_to_vertices_and_edges)
1856        psyco.bind(Domain.compute_fluxes)
1857
1858if __name__ == "__main__":
1859    pass
1860
1861
Note: See TracBrowser for help on using the repository browser.