source: branches/numpy/anuga/shallow_water/shallow_water_domain.py @ 6304

Last change on this file since 6304 was 6304, checked in by rwilson, 15 years ago

Initial commit of numpy changes. Still a long way to go.

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