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

Last change on this file since 7861 was 7861, checked in by hudson, 12 years ago

Fixed AABB parameter passing bug.

File size: 17.3 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\
185                            (domain, waveform)
186
187    Underlying domain must be specified when boundary is instantiated
188    """
189
190    def __init__(self, domain=None, function=None):
191        """ Instantiate a
192            Transmissive_n_momentum_zero_t_momentum_set_stage_boundary.
193            domain is the domain containing the boundary
194            function is the function to apply
195        """
196
197        Boundary.__init__(self)
198
199        if domain is None:
200            msg = 'Domain must be specified for this type boundary'
201            raise Exception, msg
202
203        if function is None:
204            msg = 'Function must be specified for this type boundary'
205            raise Exception, msg
206
207        self.domain = domain
208        self.function = function
209
210
211    def __repr__(self):
212        """ Return a representation of this instance. """
213        msg = 'Transmissive_n_momentum_zero_t_momentum_set_stage_boundary'
214        msg += '(%s)' % self.domain
215        return msg
216
217
218    def evaluate(self, vol_id, edge_id):
219        """Transmissive_n_momentum_zero_t_momentum_set_stage_boundary
220        return the edge momentum values of the volume they serve.
221        """
222
223        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
224
225        normal = self.domain.get_normal(vol_id, edge_id)
226
227
228        t = self.domain.get_time()
229
230        if hasattr(self.function, 'time'):
231            # Roll boundary over if time exceeds
232            while t > self.function.time[-1]:
233                msg = 'WARNING: domain time %.2f has exceeded' % t
234                msg += 'time provided in '
235                msg += 'transmissive_momentum_set_stage_boundary object.\n'
236                msg += 'I will continue, reusing the object from t==0'
237                log.critical(msg)
238                t -= self.function.time[-1]
239
240        value = self.function(t)
241        try:
242            x = float(value)
243        except:
244            x = float(value[0])
245
246        q[0] = x
247        ndotq = (normal[0]*q[1] + normal[1]*q[2])
248        q[1] = normal[0]*ndotq
249        q[2] = normal[1]*ndotq
250
251        return q
252
253class Transmissive_stage_zero_momentum_boundary(Boundary):
254    """Return same stage as those present in its neighbour volume.
255    Set momentum to zero.
256
257    Underlying domain must be specified when boundary is instantiated
258    """
259
260    def __init__(self, domain=None):
261        """ Instantiate a Transmissive (zero momentum) boundary. """
262
263        Boundary.__init__(self)
264
265        if domain is None:
266            msg = ('Domain must be specified for '
267                   'Transmissive_stage_zero_momentum boundary')
268            raise Exception, msg
269
270        self.domain = domain
271
272    def __repr__(self):
273        """ Return a representation of this instance. """
274        return 'Transmissive_stage_zero_momentum_boundary(%s)' % self.domain
275
276
277    def evaluate(self, vol_id, edge_id):
278        """Calculate transmissive (zero momentum) results. """
279
280        q = self.domain.get_conserved_quantities(vol_id, edge=edge_id)
281
282        q[1] = q[2] = 0.0
283        return q
284
285class Dirichlet_discharge_boundary(Boundary):
286    """ Class for a Dirichlet discharge boundary.
287
288    Sets stage (stage0)
289    Sets momentum (wh0) in the inward normal direction.
290
291    Underlying domain must be specified when boundary is instantiated
292    """
293
294    def __init__(self, domain=None, stage0=None, wh0=None):
295        Boundary.__init__(self)
296        """ Instantiate a Dirichlet discharge boundary.
297            domain underlying domain
298            stage0 stag
299            wh0 momentum in the inward normal direction.
300            """
301
302        if domain is None:
303            msg = 'Domain must be specified for this type of boundary'
304            raise Exception, msg
305
306        if stage0 is None:
307            raise Exception, 'Stage must be specified for this type of boundary'
308
309        if wh0 is None:
310            wh0 = 0.0
311
312        self.domain = domain
313        self.stage0 = stage0
314        self.wh0 = wh0
315
316
317    def __repr__(self):
318        """ Return a representation of this instance. """
319        return 'Dirichlet_Discharge_boundary(%s)' % self.domain
320
321    def evaluate(self, vol_id, edge_id):
322        """Set discharge in the (inward) normal direction"""
323
324        normal = self.domain.get_normal(vol_id,edge_id)
325        q = [self.stage0, -self.wh0*normal[0], -self.wh0*normal[1]]
326        return q
327
328        # FIXME: Consider this (taken from File_boundary) to allow
329        # spatial variation
330        # if vol_id is not None and edge_id is not None:
331        #     i = self.boundary_indices[ vol_id, edge_id ]
332        #     return self.F(t, point_id = i)
333        # else:
334        #     return self.F(t)
335
336
337class Inflow_boundary(Boundary):
338    """Apply given flow in m^3/s to boundary segment.
339    Depth and momentum is derived using Manning's formula.
340
341    Underlying domain must be specified when boundary is instantiated
342    """
343   
344    # FIXME (Ole): This is work in progress and definitely not finished.
345    # The associated test has been disabled
346
347    def __init__(self, domain=None, rate=0.0):
348        Boundary.__init__(self)
349
350        if domain is None:
351            msg = 'Domain must be specified for '
352            msg += 'Inflow boundary'
353            raise Exception, msg
354
355        self.domain = domain
356       
357        # FIXME(Ole): Allow rate to be time dependent as well
358        self.rate = rate
359        self.tag = None # Placeholder for tag associated with this object.
360
361    def __repr__(self):
362        return 'Inflow_boundary(%s)' %self.domain
363
364    def evaluate(self, vol_id, edge_id):
365        """Apply inflow rate at each edge of this boundary
366        """
367       
368        # First find all segments having the same tag is vol_id, edge_id
369        # This will be done the first time evaluate is called.
370        if self.tag is None:
371            boundary = self.domain.boundary
372            self.tag = boundary[(vol_id, edge_id)]       
373           
374            # Find total length of boundary with this tag
375            length = 0.0
376            for v_id, e_id in boundary:
377                if self.tag == boundary[(v_id, e_id)]:
378                    length += self.domain.mesh.get_edgelength(v_id, e_id)           
379
380            self.length = length
381            self.average_momentum = self.rate/length
382           
383           
384        # Average momentum has now been established across this boundary
385        # Compute momentum in the inward normal direction
386       
387        inward_normal = -self.domain.mesh.get_normal(vol_id, edge_id)       
388        xmomentum, ymomentum = self.average_momentum * inward_normal
389           
390        # Compute depth based on Manning's formula v = 1/n h^{2/3} sqrt(S)
391        # Where v is velocity, n is manning's coefficient, h is depth
392        # and S is the slope into the domain.
393        # Let mu be the momentum (vh), then this equation becomes:
394        #            mu = 1/n h^{5/3} sqrt(S)
395        # from which we can isolate depth to get
396        #             h = (mu n/sqrt(S) )^{3/5}
397       
398        slope = 0 # get gradient for this triangle dot normal
399       
400        # get manning coef from this triangle
401        friction = self.domain.get_quantity('friction').get_values(\
402                    location='edges', indices=[vol_id])[0]
403        mannings_n = friction[edge_id]
404
405        if slope > epsilon and mannings_n > epsilon:
406            depth = pow(self.average_momentum * mannings_n/math.sqrt(slope), \
407                        3.0/5) 
408        else:
409            depth = 1.0
410           
411        # Elevation on this edge   
412       
413        z = self.domain.get_quantity('elevation').get_values(\
414                    location='edges', indices=[vol_id])[0]
415        elevation = z[edge_id]
416           
417        # Assign conserved quantities and return
418        q = num.array([elevation + depth, xmomentum, ymomentum], num.Float)
419        return q
420
421
422       
423   
424           
425       
426class Field_boundary(Boundary):
427    """Set boundary from given field represented in an sww file containing
428    values for stage, xmomentum and ymomentum.
429
430    Optionally, the user can specify mean_stage to offset the stage provided
431    in the sww file.
432
433    This function is a thin wrapper around the generic File_boundary. The
434    difference between the file_boundary and field_boundary is only that the
435    field_boundary will allow you to change the level of the stage height when
436    you read in the boundary condition. This is very useful when running
437    different tide heights in the same area as you need only to convert one
438    boundary condition to a SWW file, ideally for tide height of 0 m
439    (saving disk space). Then you can use field_boundary to read this SWW file
440    and change the stage height (tide) on the fly depending on the scenario.
441    """
442
443    def __init__(self,
444                 filename,
445                 domain,
446                 mean_stage=0.0,
447                 time_thinning=1,
448                 time_limit=None,
449                 boundary_polygon=None,
450                 default_boundary=None,
451                 use_cache=False,
452                 verbose=False):
453        """Constructor
454
455        filename: Name of sww file containing stage and x/ymomentum
456        domain: pointer to shallow water domain for which the boundary applies
457        mean_stage: The mean water level which will be added to stage derived
458                    from the boundary condition
459        time_thinning: Will set how many time steps from the sww file read in
460                       will be interpolated to the boundary. For example if
461                       the sww file has 1 second time steps and is 24 hours
462                       in length it has 86400 time steps. If you set
463                       time_thinning to 1 it will read all these steps.
464                       If you set it to 100 it will read every 100th step eg
465                       only 864 step. This parameter is very useful to increase
466                       the speed of a model run that you are setting up
467                       and testing.
468
469        default_boundary: Must be either None or an instance of a
470                          class descending from class Boundary.
471                          This will be used in case model time exceeds
472                          that available in the underlying data.
473
474                          Note that mean_stage will also be added to this.
475
476        time_limit:
477        boundary_polygon:
478        use_cache:        True if caching is to be used.
479        verbose:          True if this method is to be verbose.
480
481        """
482
483        # Create generic file_boundary object
484        self.file_boundary = File_boundary(filename,
485                                           domain,
486                                           time_thinning=time_thinning,
487                                           time_limit=time_limit,
488                                           boundary_polygon=boundary_polygon,
489                                           default_boundary=default_boundary,
490                                           use_cache=use_cache,
491                                           verbose=verbose)
492
493        # Record information from File_boundary
494        self.F = self.file_boundary.F
495        self.domain = self.file_boundary.domain
496
497        # Record mean stage
498        self.mean_stage = mean_stage
499
500
501    def __repr__(self):
502        """ Generate a string representation of this instance. """
503        return 'Field boundary'
504
505
506    def evaluate(self, vol_id=None, edge_id=None):
507        """ Calculate 'field' boundary results.
508            vol_id and edge_id are ignored
509
510            Return linearly interpolated values based on domain.time
511        """
512
513        # Evaluate file boundary
514        q = self.file_boundary.evaluate(vol_id, edge_id)
515
516        # Adjust stage
517        for j, name in enumerate(self.domain.conserved_quantities):
518            if name == 'stage':
519                q[j] += self.mean_stage
520        return q
Note: See TracBrowser for help on using the repository browser.