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

Last change on this file since 6002 was 6002, checked in by ole, 15 years ago

Added test for Rainfall and thus also General forcing to make sure
it works with geo_referencing.

  • Property svn:keywords set to LastChangedDate LastChangedRevision LastChangedBy HeadURL Id
File size: 69.1 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: 6002 $)
70ModifiedBy:
71    $Author: ole $
72    $Date: 2008-11-26 04:39:12 +0000 (Wed, 26 Nov 2008) $
73
74"""
75
76# Subversion keywords:
77#
78# $LastChangedDate: 2008-11-26 04:39:12 +0000 (Wed, 26 Nov 2008) $
79# $LastChangedRevision: 6002 $
80# $LastChangedBy: ole $
81
82from Numeric import zeros, ones, Float, array, sum, size
83from Numeric import compress, arange
84
85from anuga.abstract_2d_finite_volumes.neighbour_mesh import segment_midpoints
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.geospatial_data.geospatial_data import ensure_geospatial
100
101from anuga.config import minimum_storable_height
102from anuga.config import minimum_allowed_height, maximum_allowed_speed
103from anuga.config import g, epsilon, beta_w, beta_w_dry,\
104     beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, tight_slope_limiters
105from anuga.config import alpha_balance
106from anuga.config import optimise_dry_cells
107from anuga.config import optimised_gradient_limiter
108from anuga.config import use_edge_limiter
109from anuga.config import use_centroid_velocities
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# Shallow water domain
121#---------------------
122class Domain(Generic_Domain):
123
124    conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
125    other_quantities = ['elevation', 'friction']
126   
127    def __init__(self,
128                 coordinates=None,
129                 vertices=None,
130                 boundary=None,
131                 tagged_elements=None,
132                 geo_reference=None,
133                 use_inscribed_circle=False,
134                 mesh_filename=None,
135                 use_cache=False,
136                 verbose=False,
137                 full_send_dict=None,
138                 ghost_recv_dict=None,
139                 processor=0,
140                 numproc=1,
141                 number_of_full_nodes=None,
142                 number_of_full_triangles=None):
143
144
145        other_quantities = ['elevation', 'friction']
146        Generic_Domain.__init__(self,
147                                coordinates,
148                                vertices,
149                                boundary,
150                                Domain.conserved_quantities,
151                                Domain.other_quantities,
152                                tagged_elements,
153                                geo_reference,
154                                use_inscribed_circle,
155                                mesh_filename,
156                                use_cache,
157                                verbose,
158                                full_send_dict,
159                                ghost_recv_dict,
160                                processor,
161                                numproc,
162                                number_of_full_nodes=number_of_full_nodes,
163                                number_of_full_triangles=number_of_full_triangles) 
164
165
166        self.set_minimum_allowed_height(minimum_allowed_height)
167       
168        self.maximum_allowed_speed = maximum_allowed_speed
169        self.g = g
170        self.beta_w      = beta_w
171        self.beta_w_dry  = beta_w_dry
172        self.beta_uh     = beta_uh
173        self.beta_uh_dry = beta_uh_dry
174        self.beta_vh     = beta_vh
175        self.beta_vh_dry = beta_vh_dry
176        self.alpha_balance = alpha_balance
177
178        self.tight_slope_limiters = tight_slope_limiters
179        self.optimise_dry_cells = optimise_dry_cells
180
181        self.forcing_terms.append(manning_friction_implicit)
182        self.forcing_terms.append(gravity)
183
184        # Stored output
185        self.store = True
186        self.format = 'sww'
187        self.set_store_vertices_uniquely(False)
188        self.minimum_storable_height = minimum_storable_height
189        self.quantities_to_be_stored = ['stage','xmomentum','ymomentum']
190
191        # Limiters
192        self.use_edge_limiter = use_edge_limiter
193        self.optimised_gradient_limiter = optimised_gradient_limiter
194        self.use_centroid_velocities = use_centroid_velocities
195
196
197    def set_all_limiters(self, beta):
198        """Shorthand to assign one constant value [0,1[ to all limiters.
199        0 Corresponds to first order, where as larger values make use of
200        the second order scheme.
201        """
202
203        self.beta_w      = beta
204        self.beta_w_dry  = beta
205        self.quantities['stage'].beta = beta
206       
207        self.beta_uh     = beta
208        self.beta_uh_dry = beta
209        self.quantities['xmomentum'].beta = beta
210       
211        self.beta_vh     = beta
212        self.beta_vh_dry = beta
213        self.quantities['ymomentum'].beta = beta
214       
215       
216
217    def set_store_vertices_uniquely(self, flag, reduction=None):
218        """Decide whether vertex values should be stored uniquely as
219        computed in the model or whether they should be reduced to one
220        value per vertex using self.reduction.
221        """
222
223        # FIXME (Ole): how about using the word continuous vertex values?
224        self.smooth = not flag
225
226        # Reduction operation for get_vertex_values
227        if reduction is None:
228            self.reduction = mean
229            #self.reduction = min  #Looks better near steep slopes
230
231
232    def set_minimum_storable_height(self, minimum_storable_height):
233        """
234        Set the minimum depth that will be recognised when writing
235        to an sww file. This is useful for removing thin water layers
236        that seems to be caused by friction creep.
237
238        The minimum allowed sww depth is in meters.
239        """
240        self.minimum_storable_height = minimum_storable_height
241
242
243    def set_minimum_allowed_height(self, minimum_allowed_height):
244        """
245        Set the minimum depth that will be recognised in the numerical
246        scheme
247
248        The minimum allowed depth is in meters.
249
250        The parameter H0 (Minimal height for flux computation)
251        is also set by this function
252        """
253
254        #FIXME (Ole): rename H0 to minimum_allowed_height_in_flux_computation
255
256        #FIXME (Ole): Maybe use histogram to identify isolated extreme speeds
257        #and deal with them adaptively similarly to how we used to use 1 order
258        #steps to recover.
259        self.minimum_allowed_height = minimum_allowed_height
260        self.H0 = minimum_allowed_height   
261       
262
263    def set_maximum_allowed_speed(self, maximum_allowed_speed):
264        """
265        Set the maximum particle speed that is allowed in water
266        shallower than minimum_allowed_height. This is useful for
267        controlling speeds in very thin layers of water and at the same time
268        allow some movement avoiding pooling of water.
269
270        """
271        self.maximum_allowed_speed = maximum_allowed_speed
272
273
274    def set_points_file_block_line_size(self,points_file_block_line_size):
275        """
276        Set the minimum depth that will be recognised when writing
277        to an sww file. This is useful for removing thin water layers
278        that seems to be caused by friction creep.
279
280        The minimum allowed sww depth is in meters.
281        """
282        self.points_file_block_line_size = points_file_block_line_size
283       
284       
285    def set_quantities_to_be_stored(self, q):
286        """Specify which quantities will be stored in the sww file.
287
288        q must be either:
289          - the name of a quantity
290          - a list of quantity names
291          - None
292
293        In the two first cases, the named quantities will be stored at
294        each yieldstep (This is in addition to the quantities elevation
295        and friction)
296       
297        If q is None, storage will be switched off altogether.
298        """
299
300
301        if q is None:
302            self.quantities_to_be_stored = []
303            self.store = False
304            return
305
306        if isinstance(q, basestring):
307            q = [q] # Turn argument into a list
308
309        # Check correcness
310        for quantity_name in q:
311            msg = 'Quantity %s is not a valid conserved quantity'\
312                  %quantity_name
313           
314            assert quantity_name in self.conserved_quantities, msg
315
316        self.quantities_to_be_stored = q
317
318
319
320    def get_wet_elements(self, indices=None):
321        """Return indices for elements where h > minimum_allowed_height
322
323        Optional argument:
324            indices is the set of element ids that the operation applies to.
325
326        Usage:
327            indices = get_wet_elements()
328
329        Note, centroid values are used for this operation           
330        """
331
332        # Water depth below which it is considered to be 0 in the model
333        # FIXME (Ole): Allow this to be specified as a keyword argument as well
334        from anuga.config import minimum_allowed_height
335
336       
337        elevation = self.get_quantity('elevation').\
338                    get_values(location='centroids', indices=indices)
339        stage = self.get_quantity('stage').\
340                get_values(location='centroids', indices=indices)       
341        depth = stage - elevation
342
343        # Select indices for which depth > 0
344        wet_indices = compress(depth > minimum_allowed_height,
345                               arange(len(depth)))
346        return wet_indices
347
348
349    def get_maximum_inundation_elevation(self, indices=None):
350        """Return highest elevation where h > 0
351
352        Optional argument:
353            indices is the set of element ids that the operation applies to.
354
355        Usage:
356            q = get_maximum_inundation_elevation()
357
358        Note, centroid values are used for this operation           
359        """
360
361        wet_elements = self.get_wet_elements(indices)
362        return self.get_quantity('elevation').\
363               get_maximum_value(indices=wet_elements)
364
365
366    def get_maximum_inundation_location(self, indices=None):
367        """Return location of highest elevation where h > 0
368
369        Optional argument:
370            indices is the set of element ids that the operation applies to.
371
372        Usage:
373            q = get_maximum_inundation_location()
374
375        Note, centroid values are used for this operation           
376        """
377
378        wet_elements = self.get_wet_elements(indices)
379        return self.get_quantity('elevation').\
380               get_maximum_location(indices=wet_elements)   
381               
382               
383               
384               
385    def get_flow_through_cross_section(self, polyline,
386                                       verbose=False):               
387        """Get the total flow through an arbitrary poly line.       
388       
389        This is a run-time equivalent of the function with same name
390        in data_manager.py
391       
392        Input:
393            polyline: Representation of desired cross section - it may contain
394                      multiple sections allowing for complex shapes. Assume
395                      absolute UTM coordinates.
396                      Format [[x0, y0], [x1, y1], ...]       
397                 
398        Output:       
399            Q: Total flow [m^3/s] across given segments.
400       
401         
402        """       
403       
404        # Find all intersections and associated triangles.
405        segments = self.get_intersecting_segments(polyline, 
406                                                  use_cache=True,
407                                                  verbose=verbose)
408
409        # Get midpoints
410        midpoints = segment_midpoints(segments)       
411       
412        # Make midpoints Geospatial instances
413        midpoints = ensure_geospatial(midpoints, self.geo_reference)       
414       
415        # Compute flow       
416        if verbose: print 'Computing flow through specified cross section'
417       
418        # Get interpolated values
419        xmomentum = self.get_quantity('xmomentum')
420        ymomentum = self.get_quantity('ymomentum')       
421       
422        uh = xmomentum.get_values(interpolation_points=midpoints, use_cache=True)
423        vh = ymomentum.get_values(interpolation_points=midpoints, use_cache=True)
424       
425        # Compute and sum flows across each segment
426        total_flow=0
427        for i in range(len(uh)):
428           
429            # Inner product of momentum vector with segment normal [m^2/s]
430            normal = segments[i].normal
431            normal_momentum = uh[i]*normal[0] + vh[i]*normal[1] 
432               
433            # Flow across this segment [m^3/s]
434            segment_flow = normal_momentum*segments[i].length
435
436            # Accumulate
437            total_flow += segment_flow
438           
439        return total_flow
440
441       
442               
443    def get_energy_through_cross_section(self, polyline,
444                                         kind='total',
445                                         verbose=False):               
446        """Obtain average energy head [m] across specified cross section.
447
448        Inputs:
449            polyline: Representation of desired cross section - it may contain
450                      multiple sections allowing for complex shapes. Assume
451                      absolute UTM coordinates.
452                      Format [[x0, y0], [x1, y1], ...]
453            kind:     Select which energy to compute.
454                      Options are 'specific' and 'total' (default)
455
456        Output:
457            E: Average energy [m] across given segments for all stored times.
458
459        The average velocity is computed for each triangle intersected by
460        the polyline and averaged weighted by segment lengths.
461
462        The typical usage of this function would be to get average energy of
463        flow in a channel, and the polyline would then be a cross section
464        perpendicular to the flow.
465
466        #FIXME (Ole) - need name for this energy reflecting that its dimension
467        is [m].
468        """
469
470        from anuga.config import g, epsilon, velocity_protection as h0       
471                                         
472       
473        # Find all intersections and associated triangles.
474        segments = self.get_intersecting_segments(polyline, 
475                                                  use_cache=True,
476                                                  verbose=verbose)
477
478        # Get midpoints
479        midpoints = segment_midpoints(segments)       
480       
481        # Make midpoints Geospatial instances
482        midpoints = ensure_geospatial(midpoints, self.geo_reference)       
483       
484        # Compute energy       
485        if verbose: print 'Computing %s energy' %kind       
486       
487        # Get interpolated values
488        stage = self.get_quantity('stage')       
489        elevation = self.get_quantity('elevation')               
490        xmomentum = self.get_quantity('xmomentum')
491        ymomentum = self.get_quantity('ymomentum')       
492
493        w = stage.get_values(interpolation_points=midpoints, use_cache=True)
494        z = elevation.get_values(interpolation_points=midpoints, use_cache=True)       
495        uh = xmomentum.get_values(interpolation_points=midpoints, use_cache=True)
496        vh = ymomentum.get_values(interpolation_points=midpoints, use_cache=True)
497        h = w-z # Depth
498       
499        # Compute total length of polyline for use with weighted averages
500        total_line_length = 0.0
501        for segment in segments:
502            total_line_length += segment.length
503       
504        # Compute and sum flows across each segment
505        average_energy=0.0
506        for i in range(len(w)):
507           
508            # Average velocity across this segment
509            if h[i] > epsilon:
510                # Use protection against degenerate velocities
511                u = uh[i]/(h[i] + h0/h[i])
512                v = vh[i]/(h[i] + h0/h[i])
513            else:
514                u = v = 0.0
515               
516            speed_squared = u*u + v*v   
517            kinetic_energy = 0.5*speed_squared/g
518           
519            if kind == 'specific':
520                segment_energy = h[i] + kinetic_energy
521            elif kind == 'total':
522                segment_energy = w[i] + kinetic_energy               
523            else:
524                msg = 'Energy kind must be either "specific" or "total".'
525                msg += ' I got %s' %kind
526               
527
528            # Add to weighted average
529            weigth = segments[i].length/total_line_length
530            average_energy += segment_energy*weigth
531             
532           
533        return average_energy
534       
535
536                       
537
538    def check_integrity(self):
539        Generic_Domain.check_integrity(self)
540
541        #Check that we are solving the shallow water wave equation
542
543        msg = 'First conserved quantity must be "stage"'
544        assert self.conserved_quantities[0] == 'stage', msg
545        msg = 'Second conserved quantity must be "xmomentum"'
546        assert self.conserved_quantities[1] == 'xmomentum', msg
547        msg = 'Third conserved quantity must be "ymomentum"'
548        assert self.conserved_quantities[2] == 'ymomentum', msg
549
550    def extrapolate_second_order_sw(self):
551        #Call correct module function
552        #(either from this module or C-extension)
553        extrapolate_second_order_sw(self)
554
555    def compute_fluxes(self):
556        #Call correct module function
557        #(either from this module or C-extension)
558        compute_fluxes(self)
559
560    def distribute_to_vertices_and_edges(self):
561        # Call correct module function
562        if self.use_edge_limiter:
563            distribute_using_edge_limiter(self)           
564        else:
565            distribute_using_vertex_limiter(self)
566
567
568
569
570    def evolve(self,
571               yieldstep = None,
572               finaltime = None,
573               duration = None,
574               skip_initial_step = False):
575        """Specialisation of basic evolve method from parent class
576        """
577
578        # Call check integrity here rather than from user scripts
579        # self.check_integrity()
580
581        msg = 'Parameter beta_w must be in the interval [0, 2['
582        assert 0 <= self.beta_w <= 2.0, msg
583
584
585        # Initial update of vertex and edge values before any STORAGE
586        # and or visualisation
587        # This is done again in the initialisation of the Generic_Domain
588        # evolve loop but we do it here to ensure the values are ok for storage
589        self.distribute_to_vertices_and_edges()
590
591        if self.store is True and self.time == 0.0:
592            self.initialise_storage()
593            # print 'Storing results in ' + self.writer.filename
594        else:
595            pass
596            # print 'Results will not be stored.'
597            # print 'To store results set domain.store = True'
598            # FIXME: Diagnostic output should be controlled by
599            # a 'verbose' flag living in domain (or in a parent class)
600
601        # Call basic machinery from parent class
602        for t in Generic_Domain.evolve(self,
603                                       yieldstep=yieldstep,
604                                       finaltime=finaltime,
605                                       duration=duration,
606                                       skip_initial_step=skip_initial_step):
607
608            # Store model data, e.g. for subsequent visualisation
609            if self.store is True:
610                self.store_timestep(self.quantities_to_be_stored)
611
612            # FIXME: Could maybe be taken from specified list
613            # of 'store every step' quantities
614
615            # Pass control on to outer loop for more specific actions
616            yield(t)
617
618
619    def initialise_storage(self):
620        """Create and initialise self.writer object for storing data.
621        Also, save x,y and bed elevation
622        """
623
624        from anuga.shallow_water.data_manager import get_dataobject
625
626        # Initialise writer
627        self.writer = get_dataobject(self, mode = 'w')
628
629        # Store vertices and connectivity
630        self.writer.store_connectivity()
631
632
633    def store_timestep(self, name):
634        """Store named quantity and time.
635
636        Precondition:
637           self.write has been initialised
638        """
639        self.writer.store_timestep(name)
640
641       
642    def timestepping_statistics(self,
643                                track_speeds=False,
644                                triangle_id=None):       
645        """Return string with time stepping statistics for printing or logging
646
647        Optional boolean keyword track_speeds decides whether to report
648        location of smallest timestep as well as a histogram and percentile
649        report.
650        """
651
652        from Numeric import sqrt
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/(sqrt(g*Vh) + epsilon)
746            Efx = Eu/(sqrt(g*Eh) + epsilon)
747            Cfx = Cu/(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/(sqrt(g*Vh) + epsilon)
763            Efy = Ev/(sqrt(g*Eh) + epsilon)
764            Cfy = Cv/(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 = zeros(3, 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
1110class Transmissive_Momentum_Set_Stage_boundary(Boundary):
1111    """Returns same momentum conserved quantities as
1112    those present in its neighbour volume.
1113    Sets stage by specifying a function f of time which may either be a
1114    vector function or a scalar function
1115
1116    Example:
1117
1118    def waveform(t):
1119        return sea_level + normalized_amplitude/cosh(t-25)**2
1120
1121    Bts = Transmissive_Momentum_Set_Stage_boundary(domain, waveform)
1122   
1123
1124    Underlying domain must be specified when boundary is instantiated
1125    """
1126
1127    def __init__(self, domain = None, function=None):
1128        Boundary.__init__(self)
1129
1130        if domain is None:
1131            msg = 'Domain must be specified for this type boundary'
1132            raise msg
1133
1134        if function is None:
1135            msg = 'Function must be specified for this type boundary'
1136            raise msg
1137
1138        self.domain   = domain
1139        self.function = function
1140
1141    def __repr__(self):
1142        return 'Transmissive_Momentum_Set_Stage_boundary(%s)' %self.domain
1143
1144    def evaluate(self, vol_id, edge_id):
1145        """Transmissive Momentum Set Stage boundaries return the edge momentum
1146        values of the volume they serve.
1147        """
1148
1149        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
1150
1151
1152        t = self.domain.get_time()
1153
1154        if hasattr(self.function, 'time'):
1155            # Roll boundary over if time exceeds           
1156            while t > self.function.time[-1]:
1157                msg = 'WARNING: domain time %.2f has exceeded' %t
1158                msg += 'time provided in '
1159                msg += 'transmissive_momentum_set_stage boundary object.\n'
1160                msg += 'I will continue, reusing the object from t==0'
1161                print msg
1162                t -= self.function.time[-1]
1163
1164
1165        value = self.function(t)
1166        try:
1167            x = float(value)
1168        except:
1169            x = float(value[0])
1170           
1171        q[0] = x
1172        return q
1173
1174
1175        # FIXME: Consider this (taken from File_boundary) to allow
1176        # spatial variation
1177        # if vol_id is not None and edge_id is not None:
1178        #     i = self.boundary_indices[ vol_id, edge_id ]
1179        #     return self.F(t, point_id = i)
1180        # else:
1181        #     return self.F(t)
1182
1183
1184
1185class Dirichlet_Discharge_boundary(Boundary):
1186    """
1187    Sets stage (stage0)
1188    Sets momentum (wh0) in the inward normal direction.
1189
1190    Underlying domain must be specified when boundary is instantiated
1191    """
1192
1193    def __init__(self, domain = None, stage0=None, wh0=None):
1194        Boundary.__init__(self)
1195
1196        if domain is None:
1197            msg = 'Domain must be specified for this type boundary'
1198            raise msg
1199
1200        if stage0 is None:
1201            raise 'set stage'
1202
1203        if wh0 is None:
1204            wh0 = 0.0
1205
1206        self.domain   = domain
1207        self.stage0  = stage0
1208        self.wh0 = wh0
1209
1210    def __repr__(self):
1211        return 'Dirichlet_Discharge_boundary(%s)' %self.domain
1212
1213    def evaluate(self, vol_id, edge_id):
1214        """Set discharge in the (inward) normal direction
1215        """
1216
1217        normal = self.domain.get_normal(vol_id,edge_id)
1218        q = [self.stage0, -self.wh0*normal[0], -self.wh0*normal[1]]
1219        return q
1220
1221
1222        # FIXME: Consider this (taken from File_boundary) to allow
1223        # spatial variation
1224        # if vol_id is not None and edge_id is not None:
1225        #     i = self.boundary_indices[ vol_id, edge_id ]
1226        #     return self.F(t, point_id = i)
1227        # else:
1228        #     return self.F(t)
1229
1230
1231class Field_boundary(Boundary):
1232    """Set boundary from given field represented in an sww file containing values
1233    for stage, xmomentum and ymomentum.
1234    Optionally, the user can specify mean_stage to offset the stage provided in the
1235    sww file.
1236
1237    This function is a thin wrapper around the generic File_boundary. The
1238    difference between the file_boundary and field_boundary is only that the
1239    field_boundary will allow you to change the level of the stage height when
1240    you read in the boundary condition. This is very useful when running
1241    different tide heights in the same area as you need only to convert one
1242    boundary condition to a SWW file, ideally for tide height of 0 m
1243    (saving disk space). Then you can use field_boundary to read this SWW file
1244    and change the stage height (tide) on the fly depending on the scenario.
1245   
1246    """
1247
1248
1249    def __init__(self, filename, domain,
1250                 mean_stage=0.0,
1251                 time_thinning=1,
1252                 boundary_polygon=None,   
1253                 default_boundary=None,                 
1254                 use_cache=False,
1255                 verbose=False):
1256        """Constructor
1257
1258        filename: Name of sww file
1259        domain: pointer to shallow water domain for which the boundary applies
1260        mean_stage: The mean water level which will be added to stage derived
1261                    from the sww file
1262        time_thinning: Will set how many time steps from the sww file read in
1263                       will be interpolated to the boundary. For example if
1264                       the sww file has 1 second time steps and is 24 hours
1265                       in length it has 86400 time steps. If you set
1266                       time_thinning to 1 it will read all these steps.
1267                       If you set it to 100 it will read every 100th step eg
1268                       only 864 step. This parameter is very useful to increase
1269                       the speed of a model run that you are setting up
1270                       and testing.
1271                       
1272        default_boundary: Must be either None or an instance of a
1273                          class descending from class Boundary.
1274                          This will be used in case model time exceeds
1275                          that available in the underlying data.
1276                                               
1277        use_cache:
1278        verbose:
1279       
1280        """
1281
1282        # Create generic file_boundary object
1283        self.file_boundary = File_boundary(filename, domain,
1284                                           time_thinning=time_thinning,
1285                                           boundary_polygon=boundary_polygon,
1286                                           default_boundary=default_boundary,
1287                                           use_cache=use_cache,
1288                                           verbose=verbose)
1289
1290       
1291        # Record information from File_boundary
1292        self.F = self.file_boundary.F
1293        self.domain = self.file_boundary.domain
1294       
1295        # Record mean stage
1296        self.mean_stage = mean_stage
1297
1298
1299    def __repr__(self):
1300        return 'Field boundary'
1301
1302
1303    def evaluate(self, vol_id=None, edge_id=None):
1304        """Return linearly interpolated values based on domain.time
1305
1306        vol_id and edge_id are ignored
1307        """
1308       
1309        # Evaluate file boundary
1310        q = self.file_boundary.evaluate(vol_id, edge_id)
1311
1312        # Adjust stage
1313        for j, name in enumerate(self.domain.conserved_quantities):
1314            if name == 'stage':
1315                q[j] += self.mean_stage
1316        return q
1317
1318   
1319
1320#-----------------------
1321# Standard forcing terms
1322#-----------------------
1323
1324def gravity(domain):
1325    """Apply gravitational pull in the presence of bed slope
1326    Wrapper calls underlying C implementation
1327    """
1328
1329    xmom = domain.quantities['xmomentum'].explicit_update
1330    ymom = domain.quantities['ymomentum'].explicit_update
1331
1332    stage = domain.quantities['stage']
1333    elevation = domain.quantities['elevation']
1334
1335    h = stage.centroid_values - elevation.centroid_values
1336    z = elevation.vertex_values
1337
1338    x = domain.get_vertex_coordinates()
1339    g = domain.g
1340   
1341
1342    from shallow_water_ext import gravity as gravity_c
1343    gravity_c(g, h, z, x, xmom, ymom) #, 1.0e-6)
1344
1345
1346
1347def manning_friction_implicit(domain):
1348    """Apply (Manning) friction to water momentum   
1349    Wrapper for c version
1350    """
1351
1352
1353    #print 'Implicit friction'
1354
1355    xmom = domain.quantities['xmomentum']
1356    ymom = domain.quantities['ymomentum']
1357
1358    w = domain.quantities['stage'].centroid_values
1359    z = domain.quantities['elevation'].centroid_values
1360
1361    uh = xmom.centroid_values
1362    vh = ymom.centroid_values
1363    eta = domain.quantities['friction'].centroid_values
1364
1365    xmom_update = xmom.semi_implicit_update
1366    ymom_update = ymom.semi_implicit_update
1367
1368    N = len(domain)
1369    eps = domain.minimum_allowed_height
1370    g = domain.g
1371
1372    from shallow_water_ext import manning_friction as manning_friction_c
1373    manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
1374
1375
1376def manning_friction_explicit(domain):
1377    """Apply (Manning) friction to water momentum   
1378    Wrapper for c version
1379    """
1380
1381    # print 'Explicit friction'
1382
1383    xmom = domain.quantities['xmomentum']
1384    ymom = domain.quantities['ymomentum']
1385
1386    w = domain.quantities['stage'].centroid_values
1387    z = domain.quantities['elevation'].centroid_values
1388
1389    uh = xmom.centroid_values
1390    vh = ymom.centroid_values
1391    eta = domain.quantities['friction'].centroid_values
1392
1393    xmom_update = xmom.explicit_update
1394    ymom_update = ymom.explicit_update
1395
1396    N = len(domain)
1397    eps = domain.minimum_allowed_height
1398    g = domain.g
1399
1400    from shallow_water_ext import manning_friction as manning_friction_c
1401    manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
1402
1403
1404# FIXME (Ole): This was implemented for use with one of the analytical solutions (Sampson?)
1405# Is it still needed (30 Oct 2007)?
1406def linear_friction(domain):
1407    """Apply linear friction to water momentum
1408
1409    Assumes quantity: 'linear_friction' to be present
1410    """
1411
1412    from math import sqrt
1413
1414    w = domain.quantities['stage'].centroid_values
1415    z = domain.quantities['elevation'].centroid_values
1416    h = w-z
1417
1418    uh = domain.quantities['xmomentum'].centroid_values
1419    vh = domain.quantities['ymomentum'].centroid_values
1420    tau = domain.quantities['linear_friction'].centroid_values
1421
1422    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1423    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1424
1425    N = len(domain) # number_of_triangles
1426    eps = domain.minimum_allowed_height
1427    g = domain.g #Not necessary? Why was this added?
1428
1429    for k in range(N):
1430        if tau[k] >= eps:
1431            if h[k] >= eps:
1432                S = -tau[k]/h[k]
1433
1434                #Update momentum
1435                xmom_update[k] += S*uh[k]
1436                ymom_update[k] += S*vh[k]
1437
1438
1439
1440#---------------------------------
1441# Experimental auxiliary functions
1442#---------------------------------
1443def check_forcefield(f):
1444    """Check that f is either
1445    1: a callable object f(t,x,y), where x and y are vectors
1446       and that it returns an array or a list of same length
1447       as x and y
1448    2: a scalar
1449    """
1450
1451    if callable(f):
1452        N = 3
1453        x = ones(3, Float)
1454        y = ones(3, Float)
1455        try:
1456            q = f(1.0, x=x, y=y)
1457        except Exception, e:
1458            msg = 'Function %s could not be executed:\n%s' %(f, e)
1459            # FIXME: Reconsider this semantics
1460            raise msg
1461
1462        try:
1463            q = array(q).astype(Float)
1464        except:
1465            msg = 'Return value from vector function %s could ' %f
1466            msg += 'not be converted into a Numeric array of floats.\n'
1467            msg += 'Specified function should return either list or array.'
1468            raise msg
1469
1470        # Is this really what we want?
1471        msg = 'Return vector from function %s ' %f
1472        msg += 'must have same lenght as input vectors'
1473        assert len(q) == N, msg
1474
1475    else:
1476        try:
1477            f = float(f)
1478        except:
1479            msg = 'Force field %s must be either a scalar' %f
1480            msg += ' or a vector function'
1481            raise Exception(msg)
1482    return f
1483
1484
1485class Wind_stress:
1486    """Apply wind stress to water momentum in terms of
1487    wind speed [m/s] and wind direction [degrees]
1488    """
1489
1490    def __init__(self, *args, **kwargs):
1491        """Initialise windfield from wind speed s [m/s]
1492        and wind direction phi [degrees]
1493
1494        Inputs v and phi can be either scalars or Python functions, e.g.
1495
1496        W = Wind_stress(10, 178)
1497
1498        #FIXME - 'normal' degrees are assumed for now, i.e. the
1499        vector (1,0) has zero degrees.
1500        We may need to convert from 'compass' degrees later on and also
1501        map from True north to grid north.
1502
1503        Arguments can also be Python functions of t,x,y as in
1504
1505        def speed(t,x,y):
1506            ...
1507            return s
1508
1509        def angle(t,x,y):
1510            ...
1511            return phi
1512
1513        where x and y are vectors.
1514
1515        and then pass the functions in
1516
1517        W = Wind_stress(speed, angle)
1518
1519        The instantiated object W can be appended to the list of
1520        forcing_terms as in
1521
1522        Alternatively, one vector valued function for (speed, angle)
1523        can be applied, providing both quantities simultaneously.
1524        As in
1525        W = Wind_stress(F), where returns (speed, angle) for each t.
1526
1527        domain.forcing_terms.append(W)
1528        """
1529
1530        from anuga.config import rho_a, rho_w, eta_w
1531        from Numeric import array, Float
1532
1533        if len(args) == 2:
1534            s = args[0]
1535            phi = args[1]
1536        elif len(args) == 1:
1537            # Assume vector function returning (s, phi)(t,x,y)
1538            vector_function = args[0]
1539            s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
1540            phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
1541        else:
1542           # Assume info is in 2 keyword arguments
1543
1544           if len(kwargs) == 2:
1545               s = kwargs['s']
1546               phi = kwargs['phi']
1547           else:
1548               raise 'Assumes two keyword arguments: s=..., phi=....'
1549
1550        self.speed = check_forcefield(s)
1551        self.phi = check_forcefield(phi)
1552
1553        self.const = eta_w*rho_a/rho_w
1554
1555
1556    def __call__(self, domain):
1557        """Evaluate windfield based on values found in domain
1558        """
1559
1560        from math import pi, cos, sin, sqrt
1561        from Numeric import Float, ones, ArrayType
1562
1563        xmom_update = domain.quantities['xmomentum'].explicit_update
1564        ymom_update = domain.quantities['ymomentum'].explicit_update
1565
1566        N = len(domain) # number_of_triangles
1567        t = domain.time
1568
1569        if callable(self.speed):
1570            xc = domain.get_centroid_coordinates()
1571            s_vec = self.speed(t, xc[:,0], xc[:,1])
1572        else:
1573            # Assume s is a scalar
1574
1575            try:
1576                s_vec = self.speed * ones(N, Float)
1577            except:
1578                msg = 'Speed must be either callable or a scalar: %s' %self.s
1579                raise msg
1580
1581
1582        if callable(self.phi):
1583            xc = domain.get_centroid_coordinates()
1584            phi_vec = self.phi(t, xc[:,0], xc[:,1])
1585        else:
1586            # Assume phi is a scalar
1587
1588            try:
1589                phi_vec = self.phi * ones(N, Float)
1590            except:
1591                msg = 'Angle must be either callable or a scalar: %s' %self.phi
1592                raise msg
1593
1594        assign_windfield_values(xmom_update, ymom_update,
1595                                s_vec, phi_vec, self.const)
1596
1597
1598def assign_windfield_values(xmom_update, ymom_update,
1599                            s_vec, phi_vec, const):
1600    """Python version of assigning wind field to update vectors.
1601    A c version also exists (for speed)
1602    """
1603    from math import pi, cos, sin, sqrt
1604
1605    N = len(s_vec)
1606    for k in range(N):
1607        s = s_vec[k]
1608        phi = phi_vec[k]
1609
1610        # Convert to radians
1611        phi = phi*pi/180
1612
1613        # Compute velocity vector (u, v)
1614        u = s*cos(phi)
1615        v = s*sin(phi)
1616
1617        # Compute wind stress
1618        S = const * sqrt(u**2 + v**2)
1619        xmom_update[k] += S*u
1620        ymom_update[k] += S*v
1621
1622
1623
1624
1625
1626class General_forcing:
1627    """General explicit forcing term for update of quantity
1628   
1629    This is used by Inflow and Rainfall for instance
1630   
1631
1632    General_forcing(quantity_name, rate, center, radius, polygon)
1633
1634    domain:     ANUGA computational domain
1635    quantity_name: Name of quantity to update.
1636                   It must be a known conserved quantity.
1637                   
1638    rate [?/s]: Total rate of change over the specified area. 
1639                This parameter can be either a constant or a
1640                function of time. Positive values indicate increases,
1641                negative values indicate decreases.
1642                Rate can be None at initialisation but must be specified
1643                before forcing term is applied (i.e. simulation has started).
1644
1645    center [m]: Coordinates at center of flow point
1646    radius [m]: Size of circular area
1647    polygon:    Arbitrary polygon
1648    default_rate: Rate to be used if rate fails (e.g. if model time exceeds its data)
1649                  Admissible types: None, constant number or function of t
1650
1651
1652    Either center, radius or polygon can be specified but not both.
1653    If neither are specified the entire domain gets updated.
1654    All coordinates to be specified in absolute UTM coordinates (x, y) assuming the zone of domain.   
1655   
1656    Inflow or Rainfall for examples of use
1657    """
1658
1659
1660    # FIXME (AnyOne) : Add various methods to allow spatial variations
1661
1662    def __init__(self,
1663                 domain,
1664                 quantity_name,
1665                 rate=0.0,
1666                 center=None, radius=None,
1667                 polygon=None,
1668                 default_rate=None,
1669                 verbose=False):
1670                     
1671        if center is None:
1672            msg = 'I got radius but no center.'       
1673            assert radius is None, msg
1674           
1675        if radius is None:
1676            msg += 'I got center but no radius.'       
1677            assert center is None, msg
1678           
1679           
1680                     
1681        from math import pi, cos, sin
1682
1683        self.domain = domain
1684        self.quantity_name = quantity_name
1685        self.rate = rate
1686        self.center = ensure_numeric(center)
1687        self.radius = radius
1688        self.polygon = polygon       
1689        self.verbose = verbose
1690        self.value = 0.0 # Can be used to remember value at
1691                         # previous timestep in order to obtain rate
1692
1693        bounding_polygon = domain.get_boundary_polygon() # Returns absolute coordinates
1694
1695
1696        # Update area if applicable
1697        self.exchange_area = None       
1698        if center is not None and radius is not None:
1699            assert len(center) == 2
1700            msg = 'Polygon cannot be specified when center and radius are'
1701            assert polygon is None, msg
1702
1703            self.exchange_area = radius**2*pi
1704
1705            # Check that circle center lies within the mesh.
1706            msg = 'Center %s specified for forcing term did not' %(str(center))
1707            msg += 'fall within the domain boundary.'
1708            assert is_inside_polygon(center, bounding_polygon), msg
1709
1710            # Check that circle periphery lies within the mesh.
1711            N = 100
1712            periphery_points = []
1713            for i in range(N):
1714
1715                theta = 2*pi*i/100
1716               
1717                x = center[0] + radius*cos(theta)
1718                y = center[1] + radius*sin(theta)
1719
1720                periphery_points.append([x,y])
1721
1722
1723            for point in periphery_points:
1724                msg = 'Point %s on periphery for forcing term' %(str(point))
1725                msg += ' did not fall within the domain boundary.'
1726                assert is_inside_polygon(point, bounding_polygon), msg
1727
1728       
1729        if polygon is not None:
1730            self.exchange_area = polygon_area(self.polygon)
1731
1732            # Check that polygon lies within the mesh.
1733            for point in self.polygon:
1734                msg = 'Point %s in polygon for forcing term' %(str(point))
1735                msg += ' did not fall within the domain boundary.'
1736                assert is_inside_polygon(point, bounding_polygon), msg
1737
1738
1739        # Pointer to update vector
1740        self.update = domain.quantities[self.quantity_name].explicit_update
1741
1742        # Determine indices in flow area
1743        N = len(domain)   
1744        points = domain.get_centroid_coordinates(absolute=True)
1745
1746        # Calculate indices in exchange area for this forcing term
1747        self.exchange_indices = None
1748        if self.center is not None and self.radius is not None:
1749            # Inlet is circular
1750           
1751            inlet_region = 'center=%s, radius=%s' %(self.center, self.radius) 
1752           
1753            self.exchange_indices = []
1754            for k in range(N):
1755                x, y = points[k,:] # Centroid
1756               
1757                c = self.center
1758                if ((x-c[0])**2+(y-c[1])**2) < self.radius**2:
1759                    self.exchange_indices.append(k)
1760                   
1761        if self.polygon is not None:                   
1762            # Inlet is polygon
1763           
1764            inlet_region = 'polygon=%s' %(self.polygon) 
1765                       
1766            self.exchange_indices = inside_polygon(points, self.polygon)
1767           
1768           
1769        if self.exchange_indices is not None:
1770            #print inlet_region
1771       
1772            if len(self.exchange_indices) == 0:
1773                msg = 'No triangles have been identified in '
1774                msg += 'specified region: %s' %inlet_region
1775                raise Exception, msg
1776           
1777        # Check and store default_rate
1778        msg = 'Keyword argument default_rate must be either None ' 
1779        msg += 'or a function of time.\n I got %s' %(str(default_rate)) 
1780        assert default_rate is None or \
1781               type(default_rate) in [IntType, FloatType] or \
1782               callable(default_rate), msg
1783       
1784        if default_rate is not None:
1785
1786            # If it is a constant, make it a function
1787            if not callable(default_rate):
1788                tmp = default_rate
1789                default_rate = lambda t: tmp
1790
1791           
1792            # Check that default_rate is a function of one argument
1793            try:
1794                default_rate(0.0)
1795            except:
1796                raise Exception, msg
1797
1798        self.default_rate = default_rate
1799        self.default_rate_invoked = False    # Flag       
1800       
1801
1802    def __call__(self, domain):
1803        """Apply inflow function at time specified in domain and update stage
1804        """
1805
1806        # Call virtual method allowing local modifications
1807
1808        t = domain.get_time()
1809        try:
1810            rate = self.update_rate(t)
1811        except Modeltime_too_early, e:
1812            raise Modeltime_too_early, e
1813        except Modeltime_too_late, e:
1814            if self.default_rate is None:
1815                raise Exception, e # Reraise exception
1816            else:
1817                # Pass control to default rate function
1818                rate = self.default_rate(t)
1819               
1820                if self.default_rate_invoked is False:
1821                    # Issue warning the first time
1822                    msg = '%s' %str(e)
1823                    msg += 'Instead I will use the default rate: %s\n'\
1824                        %str(self.default_rate) 
1825                    msg += 'Note: Further warnings will be supressed'
1826                    warn(msg)
1827                   
1828                    # FIXME (Ole): Replace this crude flag with
1829                    # Python's ability to print warnings only once.
1830                    # See http://docs.python.org/lib/warning-filter.html
1831                    self.default_rate_invoked = True
1832                   
1833
1834           
1835               
1836
1837        if rate is None:
1838            msg = 'Attribute rate must be specified in General_forcing'
1839            msg += ' or its descendants before attempting to call it'
1840            raise Exception, msg
1841       
1842
1843        # Now rate is a number
1844        if self.verbose is True:
1845            print 'Rate of %s at time = %.2f = %f' %(self.quantity_name,
1846                                                     domain.get_time(),
1847                                                     rate)
1848
1849
1850        if self.exchange_indices is None:
1851            self.update[:] += rate
1852        else:
1853            # Brute force assignment of restricted rate
1854            for k in self.exchange_indices:
1855                self.update[k] += rate
1856
1857
1858    def update_rate(self, t):
1859        """Virtual method allowing local modifications by writing an
1860        overriding version in descendant
1861       
1862        """
1863        if callable(self.rate):
1864            rate = self.rate(t)
1865        else:
1866            rate = self.rate
1867
1868        return rate
1869
1870
1871    def get_quantity_values(self):
1872        """Return values for specified quantity restricted to opening
1873        """
1874       
1875        q = self.domain.quantities[self.quantity_name]
1876        return q.get_values(indices=self.exchange_indices)
1877   
1878
1879    def set_quantity_values(self, val):
1880        """Set values for specified quantity restricted to opening
1881        """
1882
1883        q = self.domain.quantities[self.quantity_name]               
1884        q.set_values(val, indices=self.exchange_indices)   
1885
1886
1887
1888class Rainfall(General_forcing):
1889    """Class Rainfall - general 'rain over entire domain' forcing term.
1890   
1891    Used for implementing Rainfall over the entire domain.
1892       
1893        Current Limited to only One Gauge..
1894       
1895        Need to add Spatial Varying Capability
1896        (This module came from copying and amending the Inflow Code)
1897   
1898    Rainfall(rain)
1899
1900    domain   
1901    rain [mm/s]:  Total rain rate over the specified domain. 
1902                  NOTE: Raingauge Data needs to reflect the time step.
1903                  IE: if Gauge is mm read at a time step, then the input
1904                  here is as mm/(timeStep) so 10mm in 5minutes becomes
1905                  10/(5x60) = 0.0333mm/s.
1906       
1907       
1908                  This parameter can be either a constant or a
1909                  function of time. Positive values indicate inflow,
1910                  negative values indicate outflow.
1911                  (and be used for Infiltration - Write Seperate Module)
1912                  The specified flow will be divided by the area of
1913                  the inflow region and then applied to update the
1914                  stage quantity.
1915
1916    polygon: Specifies a polygon to restrict the rainfall.
1917   
1918    Examples
1919    How to put them in a run File...
1920       
1921    #------------------------------------------------------------------------
1922    # Setup specialised forcing terms
1923    #------------------------------------------------------------------------
1924    # This is the new element implemented by Ole and Rudy to allow direct
1925    # input of Inflow in mm/s
1926
1927    catchmentrainfall = Rainfall(rain=file_function('Q100_2hr_Rain.tms')) 
1928                        # Note need path to File in String.
1929                        # Else assumed in same directory
1930
1931    domain.forcing_terms.append(catchmentrainfall)
1932    """
1933
1934   
1935    def __init__(self,
1936                 domain,
1937                 rate=0.0,
1938                 center=None, radius=None,
1939                 polygon=None,
1940                 default_rate=None,                 
1941                 verbose=False):
1942
1943        # Converting mm/s to m/s to apply in ANUGA)
1944        if callable(rate):
1945            rain = lambda t: rate(t)/1000.0
1946        else:
1947            rain = rate/1000.0
1948
1949        if default_rate is not None:   
1950            if callable(default_rate):
1951                default_rain = lambda t: default_rate(t)/1000.0
1952            else:
1953                default_rain = default_rate/1000.0
1954        else:
1955            default_rain = None
1956
1957
1958           
1959           
1960        General_forcing.__init__(self,
1961                                 domain,
1962                                 'stage',
1963                                 rate=rain,
1964                                 center=center, radius=radius,
1965                                 polygon=polygon,
1966                                 default_rate=default_rain,
1967                                 verbose=verbose)
1968
1969       
1970
1971
1972
1973
1974class Inflow(General_forcing):
1975    """Class Inflow - general 'rain and drain' forcing term.
1976   
1977    Useful for implementing flows in and out of the domain.
1978   
1979    Inflow(flow, center, radius, polygon)
1980
1981    domain
1982    rate [m^3/s]: Total flow rate over the specified area. 
1983                  This parameter can be either a constant or a
1984                  function of time. Positive values indicate inflow,
1985                  negative values indicate outflow.
1986                  The specified flow will be divided by the area of
1987                  the inflow region and then applied to update stage.     
1988    center [m]: Coordinates at center of flow point
1989    radius [m]: Size of circular area
1990    polygon:    Arbitrary polygon.
1991
1992    Either center, radius or polygon must be specified
1993   
1994    Examples
1995
1996    # Constant drain at 0.003 m^3/s.
1997    # The outflow area is 0.07**2*pi=0.0154 m^2
1998    # This corresponds to a rate of change of 0.003/0.0154 = 0.2 m/s
1999    #                                     
2000    Inflow((0.7, 0.4), 0.07, -0.003)
2001
2002
2003    # Tap turning up to a maximum inflow of 0.0142 m^3/s.
2004    # The inflow area is 0.03**2*pi = 0.00283 m^2
2005    # This corresponds to a rate of change of 0.0142/0.00283 = 5 m/s     
2006    # over the specified area
2007    Inflow((0.5, 0.5), 0.03, lambda t: min(0.01*t, 0.0142))
2008
2009
2010    #------------------------------------------------------------------------
2011    # Setup specialised forcing terms
2012    #------------------------------------------------------------------------
2013    # This is the new element implemented by Ole to allow direct input
2014    # of Inflow in m^3/s
2015
2016    hydrograph = Inflow(center=(320, 300), radius=10,
2017                        rate=file_function('Q/QPMF_Rot_Sub13.tms'))
2018
2019    domain.forcing_terms.append(hydrograph)
2020   
2021    """
2022
2023
2024    def __init__(self,
2025                 domain,
2026                 rate=0.0,
2027                 center=None, radius=None,
2028                 polygon=None,
2029                 default_rate=None,
2030                 verbose=False):                 
2031
2032
2033        # Create object first to make area is available
2034        General_forcing.__init__(self,
2035                                 domain,
2036                                 'stage',
2037                                 rate=rate,
2038                                 center=center, radius=radius,
2039                                 polygon=polygon,
2040                                 default_rate=default_rate,
2041                                 verbose=verbose)
2042
2043    def update_rate(self, t):
2044        """Virtual method allowing local modifications by writing an
2045        overriding version in descendant
2046
2047        This one converts m^3/s to m/s which can be added directly
2048        to 'stage' in ANUGA
2049        """
2050
2051        if callable(self.rate):
2052            _rate = self.rate(t)/self.exchange_area
2053        else:
2054            _rate = self.rate/self.exchange_area
2055
2056        return _rate
2057
2058
2059
2060
2061#------------------
2062# Initialise module
2063#------------------
2064
2065
2066from anuga.utilities import compile
2067if compile.can_use_C_extension('shallow_water_ext.c'):
2068    # Underlying C implementations can be accessed
2069
2070    from shallow_water_ext import rotate, assign_windfield_values
2071else:
2072    msg = 'C implementations could not be accessed by %s.\n ' %__file__
2073    msg += 'Make sure compile_all.py has been run as described in '
2074    msg += 'the ANUGA installation guide.'
2075    raise Exception, msg
2076
2077
2078# Optimisation with psyco
2079from anuga.config import use_psyco
2080if use_psyco:
2081    try:
2082        import psyco
2083    except:
2084        import os
2085        if os.name == 'posix' and os.uname()[4] in ['x86_64', 'ia64']:
2086            pass
2087            #Psyco isn't supported on 64 bit systems, but it doesn't matter
2088        else:
2089            msg = 'WARNING: psyco (speedup) could not import'+\
2090                  ', you may want to consider installing it'
2091            print msg
2092    else:
2093        psyco.bind(Domain.distribute_to_vertices_and_edges)
2094        psyco.bind(Domain.compute_fluxes)
2095
2096if __name__ == "__main__":
2097    pass
2098
2099
Note: See TracBrowser for help on using the repository browser.