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

Last change on this file since 8920 was 8920, checked in by steve, 12 years ago

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