source: trunk/anuga_core/source/anuga/shallow_water/boundaries.py @ 7765

Last change on this file since 7765 was 7765, checked in by hudson, 13 years ago

Refactoring - moved file conversion routines to file_conversion folder, moved file loading/saving functions to file module.

File size: 17.6 KB
Line 
1"""
2Boundary conditions - specific to the shallow water wave equation
3
4Title: ANUGA boundaries with dependencies on shallow_water_domain
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: 2010
12
13Description:
14    This module contains boundary functions for ANUGA that are specific
15    to the shallow water Domain class.
16   
17Constraints: See GPL license in the user guide
18Version: 1.0 ($Revision: 7731 $)
19ModifiedBy:
20    $Author: hudson $
21    $Date: 2010-05-18 14:54:05 +1000 (Tue, 18 May 2010) $
22"""
23
24
25from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
26     import Boundary, File_boundary
27import numpy as num
28     
29
30from anuga.utilities import compile
31if compile.can_use_C_extension('shallow_water_ext.c'):
32    # Underlying C implementations can be accessed
33    from shallow_water_ext import rotate
34else:
35    msg = 'C implementations could not be accessed by %s.\n ' % __file__
36    msg += 'Make sure compile_all.py has been run as described in '
37    msg += 'the ANUGA installation guide.'
38    raise Exception, msg
39
40
41class Reflective_boundary(Boundary):
42    """Reflective boundary returns same conserved quantities as
43    those present in its neighbour volume but reflected.
44
45    This class is specific to the shallow water equation as it
46    works with the momentum quantities assumed to be the second
47    and third conserved quantities.
48    """
49
50    ##
51    # @brief Instantiate a Reflective_boundary.
52    # @param domain
53    def __init__(self, domain=None):
54        Boundary.__init__(self)
55
56        if domain is None:
57            msg = 'Domain must be specified for reflective boundary'
58            raise Exception, msg
59
60        # Handy shorthands
61        self.stage = domain.quantities['stage'].edge_values
62        self.xmom = domain.quantities['xmomentum'].edge_values
63        self.ymom = domain.quantities['ymomentum'].edge_values
64        self.normals = domain.normals
65
66        self.conserved_quantities = num.zeros(3, num.float)
67
68    ##
69    # @brief Return a representation of this instance.
70    def __repr__(self):
71        return 'Reflective_boundary'
72
73    ##
74    # @brief Calculate reflections (reverse outward momentum).
75    # @param vol_id
76    # @param edge_id
77    def evaluate(self, vol_id, edge_id):
78        """Reflective boundaries reverses the outward momentum
79        of the volume they serve.
80        """
81
82        q = self.conserved_quantities
83        q[0] = self.stage[vol_id, edge_id]
84        q[1] = self.xmom[vol_id, edge_id]
85        q[2] = self.ymom[vol_id, edge_id]
86
87        normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
88
89        r = rotate(q, normal, direction = 1)
90        r[1] = -r[1]
91        q = rotate(r, normal, direction = -1)
92
93        return q
94
95
96class Transmissive_momentum_set_stage_boundary(Boundary):
97    """Returns same momentum conserved quantities as
98    those present in its neighbour volume.
99    Sets stage by specifying a function f of time which may either be a
100    vector function or a scalar function
101
102    Example:
103
104    def waveform(t):
105        return sea_level + normalized_amplitude/cosh(t-25)**2
106
107    Bts = Transmissive_momentum_set_stage_boundary(domain, waveform)
108
109    Underlying domain must be specified when boundary is instantiated
110    """
111
112    def __init__(self, domain=None, function=None):
113        Boundary.__init__(self)
114        """ Instantiate a Transmissive_momentum_set_stage_boundary. """
115
116
117        if domain is None:
118            msg = 'Domain must be specified for this type boundary'
119            raise Exception, msg
120
121        if function is None:
122            msg = 'Function must be specified for this type boundary'
123            raise Exception, msg
124
125        self.domain = domain
126        self.function = function
127
128
129    def __repr__(self):
130        """ Return a representation of this instance. """
131        return 'Transmissive_momentum_set_stage_boundary(%s)' % self.domain
132
133    def evaluate(self, vol_id, edge_id):
134        """Transmissive momentum set stage boundaries return the edge momentum
135        values of the volume they serve.
136
137        vol_id is volume id
138        edge_id is the edge within the volume
139        """
140
141        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
142        t = self.domain.get_time()
143
144        if hasattr(self.function, 'time'):
145            # Roll boundary over if time exceeds
146            while t > self.function.time[-1]:
147                msg = 'WARNING: domain time %.2f has exceeded' % t
148                msg += 'time provided in '
149                msg += 'transmissive_momentum_set_stage_boundary object.\n'
150                msg += 'I will continue, reusing the object from t==0'
151                log.critical(msg)
152                t -= self.function.time[-1]
153
154        value = self.function(t)
155        try:
156            x = float(value)
157        except:
158            x = float(value[0])
159
160        q[0] = x
161           
162        return q
163
164        # FIXME: Consider this (taken from File_boundary) to allow
165        # spatial variation
166        # if vol_id is not None and edge_id is not None:
167        #     i = self.boundary_indices[ vol_id, edge_id ]
168        #     return self.F(t, point_id = i)
169        # else:
170        #     return self.F(t)
171
172
173class Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(Boundary):
174    """Returns the same normal momentum as that
175    present in neighbour volume edge. Zero out the tangential momentum.
176    Sets stage by specifying a function f of time which may either be a
177    vector function or a scalar function
178
179    Example:
180
181    def waveform(t):
182        return sea_level + normalized_amplitude/cosh(t-25)**2
183
184    Bts = Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(domain, waveform)
185
186    Underlying domain must be specified when boundary is instantiated
187    """
188
189    def __init__(self, domain=None, function=None):
190        """ Instantiate a
191            Transmissive_n_momentum_zero_t_momentum_set_stage_boundary.
192            domain is the domain containing the boundary
193            function is the function to apply
194        """
195
196        Boundary.__init__(self)
197
198        if domain is None:
199            msg = 'Domain must be specified for this type boundary'
200            raise Exception, msg
201
202        if function is None:
203            msg = 'Function must be specified for this type boundary'
204            raise Exception, msg
205
206        self.domain = domain
207        self.function = function
208
209
210    def __repr__(self):
211        """ Return a representation of this instance. """
212        msg = 'Transmissive_n_momentum_zero_t_momentum_set_stage_boundary'
213        msg += '(%s)' % self.domain
214        return msg
215
216
217    def evaluate(self, vol_id, edge_id):
218        """Transmissive_n_momentum_zero_t_momentum_set_stage_boundary
219        return the edge momentum values of the volume they serve.
220        """
221
222        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
223
224        normal = self.domain.get_normal(vol_id, edge_id)
225
226
227        t = self.domain.get_time()
228
229        if hasattr(self.function, 'time'):
230            # Roll boundary over if time exceeds
231            while t > self.function.time[-1]:
232                msg = 'WARNING: domain time %.2f has exceeded' % t
233                msg += 'time provided in '
234                msg += 'transmissive_momentum_set_stage_boundary object.\n'
235                msg += 'I will continue, reusing the object from t==0'
236                log.critical(msg)
237                t -= self.function.time[-1]
238
239        value = self.function(t)
240        try:
241            x = float(value)
242        except:
243            x = float(value[0])
244
245        ## import math
246        ## if vol_id == 9433:
247        ##     print 'vol_id = ',vol_id, ' edge_id = ',edge_id, 'q = ', q, ' x = ',x
248        ##     print 'normal = ', normal
249        ##     print 'n . p = ', (normal[0]*q[1] + normal[1]*q[2])
250        ##     print 't . p = ', (normal[1]*q[1] - normal[0]*q[2])
251
252
253        q[0] = x
254        ndotq = (normal[0]*q[1] + normal[1]*q[2])
255        q[1] = normal[0]*ndotq
256        q[2] = normal[1]*ndotq
257
258        return q
259
260class Transmissive_stage_zero_momentum_boundary(Boundary):
261    """Return same stage as those present in its neighbour volume.
262    Set momentum to zero.
263
264    Underlying domain must be specified when boundary is instantiated
265    """
266
267    def __init__(self, domain=None):
268        """ Instantiate a Transmissive (zero momentum) boundary. """
269
270        Boundary.__init__(self)
271
272        if domain is None:
273            msg = ('Domain must be specified for '
274                   'Transmissive_stage_zero_momentum boundary')
275            raise Exception, msg
276
277        self.domain = domain
278
279    def __repr__(self):
280        """ Return a representation of this instance. """
281        return 'Transmissive_stage_zero_momentum_boundary(%s)' % self.domain
282
283
284    def evaluate(self, vol_id, edge_id):
285        """Calculate transmissive (zero momentum) results. """
286
287        q = self.domain.get_conserved_quantities(vol_id, edge=edge_id)
288
289        q[1] = q[2] = 0.0
290        return q
291
292class Dirichlet_discharge_boundary(Boundary):
293    """ Class for a Dirichlet discharge boundary.
294
295    Sets stage (stage0)
296    Sets momentum (wh0) in the inward normal direction.
297
298    Underlying domain must be specified when boundary is instantiated
299    """
300
301    def __init__(self, domain=None, stage0=None, wh0=None):
302        Boundary.__init__(self)
303        """ Instantiate a Dirichlet discharge boundary.
304            domain underlying domain
305            stage0 stag
306            wh0 momentum in the inward normal direction.
307            """
308
309        if domain is None:
310            msg = 'Domain must be specified for this type of boundary'
311            raise Exception, msg
312
313        if stage0 is None:
314            raise Exception, 'Stage must be specified for this type of boundary'
315
316        if wh0 is None:
317            wh0 = 0.0
318
319        self.domain = domain
320        self.stage0 = stage0
321        self.wh0 = wh0
322
323
324    def __repr__(self):
325        """ Return a representation of this instance. """
326        return 'Dirichlet_Discharge_boundary(%s)' % self.domain
327
328    def evaluate(self, vol_id, edge_id):
329        """Set discharge in the (inward) normal direction"""
330
331        normal = self.domain.get_normal(vol_id,edge_id)
332        q = [self.stage0, -self.wh0*normal[0], -self.wh0*normal[1]]
333        return q
334
335        # FIXME: Consider this (taken from File_boundary) to allow
336        # spatial variation
337        # if vol_id is not None and edge_id is not None:
338        #     i = self.boundary_indices[ vol_id, edge_id ]
339        #     return self.F(t, point_id = i)
340        # else:
341        #     return self.F(t)
342
343
344class Inflow_boundary(Boundary):
345    """Apply given flow in m^3/s to boundary segment.
346    Depth and momentum is derived using Manning's formula.
347
348    Underlying domain must be specified when boundary is instantiated
349    """
350   
351    # FIXME (Ole): This is work in progress and definitely not finished.
352    # The associated test has been disabled
353
354    def __init__(self, domain=None, rate=0.0):
355        Boundary.__init__(self)
356
357        if domain is None:
358            msg = 'Domain must be specified for '
359            msg += 'Inflow boundary'
360            raise Exception, msg
361
362        self.domain = domain
363       
364        # FIXME(Ole): Allow rate to be time dependent as well
365        self.rate = rate
366        self.tag = None # Placeholder for tag associated with this object.
367
368    def __repr__(self):
369        return 'Inflow_boundary(%s)' %self.domain
370
371    def evaluate(self, vol_id, edge_id):
372        """Apply inflow rate at each edge of this boundary
373        """
374       
375        # First find all segments having the same tag is vol_id, edge_id
376        # This will be done the first time evaluate is called.
377        if self.tag is None:
378            boundary = self.domain.boundary
379            self.tag = boundary[(vol_id, edge_id)]       
380           
381            # Find total length of boundary with this tag
382            length = 0.0
383            for v_id, e_id in boundary:
384                if self.tag == boundary[(v_id, e_id)]:
385                    length += self.domain.mesh.get_edgelength(v_id, e_id)           
386
387            self.length = length
388            self.average_momentum = self.rate/length
389           
390           
391        # Average momentum has now been established across this boundary
392        # Compute momentum in the inward normal direction
393       
394        inward_normal = -self.domain.mesh.get_normal(vol_id, edge_id)       
395        xmomentum, ymomentum = self.average_momentum * inward_normal
396           
397        # Compute depth based on Manning's formula v = 1/n h^{2/3} sqrt(S)
398        # Where v is velocity, n is manning's coefficient, h is depth
399        # and S is the slope into the domain.
400        # Let mu be the momentum (vh), then this equation becomes: mu = 1/n h^{5/3} sqrt(S)
401        # from which we can isolate depth to get
402        # h = (mu n/sqrt(S) )^{3/5}
403       
404        slope = 0 # get gradient for this triangle dot normal
405       
406        # get manning coef from this triangle
407        friction = self.domain.get_quantity('friction').get_values(location='edges', 
408                                                                   indices=[vol_id])[0]
409        mannings_n = friction[edge_id]
410
411        if slope > epsilon and mannings_n > epsilon:
412            depth = pow(self.average_momentum * mannings_n/math.sqrt(slope), 3.0/5) 
413        else:
414            depth = 1.0
415           
416        # Elevation on this edge   
417       
418        z = self.domain.get_quantity('elevation').get_values(location='edges', 
419                                                             indices=[vol_id])[0]
420        elevation = z[edge_id]
421           
422        # Assign conserved quantities and return
423        q = num.array([elevation + depth, xmomentum, ymomentum], num.Float)
424        return q
425
426
427       
428   
429           
430       
431class Field_boundary(Boundary):
432    """Set boundary from given field represented in an sww file containing
433    values for stage, xmomentum and ymomentum.
434
435    Optionally, the user can specify mean_stage to offset the stage provided
436    in the sww file.
437
438    This function is a thin wrapper around the generic File_boundary. The
439    difference between the file_boundary and field_boundary is only that the
440    field_boundary will allow you to change the level of the stage height when
441    you read in the boundary condition. This is very useful when running
442    different tide heights in the same area as you need only to convert one
443    boundary condition to a SWW file, ideally for tide height of 0 m
444    (saving disk space). Then you can use field_boundary to read this SWW file
445    and change the stage height (tide) on the fly depending on the scenario.
446    """
447
448    def __init__(self,
449                 filename,
450                 domain,
451                 mean_stage=0.0,
452                 time_thinning=1,
453                 time_limit=None,
454                 boundary_polygon=None,
455                 default_boundary=None,
456                 use_cache=False,
457                 verbose=False):
458        """Constructor
459
460        filename: Name of sww file containing stage and x/ymomentum
461        domain: pointer to shallow water domain for which the boundary applies
462        mean_stage: The mean water level which will be added to stage derived
463                    from the boundary condition
464        time_thinning: Will set how many time steps from the sww file read in
465                       will be interpolated to the boundary. For example if
466                       the sww file has 1 second time steps and is 24 hours
467                       in length it has 86400 time steps. If you set
468                       time_thinning to 1 it will read all these steps.
469                       If you set it to 100 it will read every 100th step eg
470                       only 864 step. This parameter is very useful to increase
471                       the speed of a model run that you are setting up
472                       and testing.
473
474        default_boundary: Must be either None or an instance of a
475                          class descending from class Boundary.
476                          This will be used in case model time exceeds
477                          that available in the underlying data.
478
479                          Note that mean_stage will also be added to this.
480
481        time_limit:
482        boundary_polygon:
483        use_cache:        True if caching is to be used.
484        verbose:          True if this method is to be verbose.
485
486        """
487
488        # Create generic file_boundary object
489        self.file_boundary = File_boundary(filename,
490                                           domain,
491                                           time_thinning=time_thinning,
492                                           time_limit=time_limit,
493                                           boundary_polygon=boundary_polygon,
494                                           default_boundary=default_boundary,
495                                           use_cache=use_cache,
496                                           verbose=verbose)
497
498        # Record information from File_boundary
499        self.F = self.file_boundary.F
500        self.domain = self.file_boundary.domain
501
502        # Record mean stage
503        self.mean_stage = mean_stage
504
505
506    def __repr__(self):
507        """ Generate a string representation of this instance. """
508        return 'Field boundary'
509
510
511    def evaluate(self, vol_id=None, edge_id=None):
512        """ Calculate 'field' boundary results.
513            vol_id and edge_id are ignored
514
515            Return linearly interpolated values based on domain.time
516        """
517
518        # Evaluate file boundary
519        q = self.file_boundary.evaluate(vol_id, edge_id)
520
521        # Adjust stage
522        for j, name in enumerate(self.domain.conserved_quantities):
523            if name == 'stage':
524                q[j] += self.mean_stage
525        return q
Note: See TracBrowser for help on using the repository browser.