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

Last change on this file since 9558 was 9374, checked in by davies, 10 years ago

Sww merge using even less memory + slight boundary condition cleanup
File size: 32.7 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        Gareth Davies, gareth.davies.ga.code@gmail.com
11
12CreationDate: 2010
13
14Description:
15    This module contains boundary functions for ANUGA that are specific
16    to the shallow water Domain class.
17   
18Constraints: See GPL license in the user guide
19Version: 1.0 ($Revision: 7731 $)
20ModifiedBy:
21    $Author: hudson $
22    $Date: 2010-05-18 14:54:05 +1000 (Tue, 18 May 2010) $
23"""
24
25
26from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
27     import Boundary, File_boundary
28import numpy as num
29
30import anuga.utilities.log as log
31from anuga.fit_interpolate.interpolate import Modeltime_too_late
32from anuga.fit_interpolate.interpolate import Modeltime_too_early
33from anuga.config import g as gravity
34     
35from shallow_water_ext import rotate
36
37#from anuga.utilities import compile
38#if compile.can_use_C_extension('shallow_water_ext.c'):
39#    # Underlying C implementations can be accessed
40#    from shallow_water_ext import rotate
41#else:
42#    msg = 'C implementations could not be accessed by %s.\n ' % __file__
43#    msg += 'Make sure compile_all.py has been run as described in '
44#    msg += 'the ANUGA installation guide.'
45#    raise Exception, msg
46
47
48class Reflective_boundary(Boundary):
49    """Reflective boundary returns same conserved quantities as
50    those present in its neighbour volume but reflected.
51
52    This class is specific to the shallow water equation as it
53    works with the momentum quantities assumed to be the second
54    and third conserved quantities.
55    """
56
57    def __init__(self, domain=None):
58        Boundary.__init__(self)
59
60        if domain is None:
61            msg = 'Domain must be specified for reflective boundary'
62            raise Exception, msg
63
64        # Handy shorthands
65        self.stage = domain.quantities['stage'].edge_values
66        self.xmom = domain.quantities['xmomentum'].edge_values
67        self.ymom = domain.quantities['ymomentum'].edge_values
68        self.normals = domain.normals
69
70        self.conserved_quantities = num.zeros(3, num.float)
71
72    def __repr__(self):
73        return 'Reflective_boundary'
74
75    def evaluate(self, vol_id, edge_id):
76        """Calculate reflections (reverse outward momentum).
77
78        vol_id   
79        edge_id 
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
96
97    def evaluate_segment(self, domain, segment_edges):
98        """Apply reflective BC on the boundary edges defined by
99        segment_edges
100        """
101
102        if segment_edges is None:
103            return
104        if domain is None:
105            return
106
107
108        ids = segment_edges
109        vol_ids  = domain.boundary_cells[ids]
110        edge_ids = domain.boundary_edges[ids]
111
112        Stage = domain.quantities['stage']
113        Elev  = domain.quantities['elevation']
114        Height= domain.quantities['height']
115        Xmom  = domain.quantities['xmomentum']
116        Ymom  = domain.quantities['ymomentum']
117        Xvel  = domain.quantities['xvelocity']
118        Yvel  = domain.quantities['yvelocity']
119
120        Normals = domain.normals
121
122        #print vol_ids
123        #print edge_ids
124        #Normals.reshape((4,3,2))
125        #print Normals.shape
126        #print Normals[vol_ids, 2*edge_ids]
127        #print Normals[vol_ids, 2*edge_ids+1]
128       
129        n1  = Normals[vol_ids,2*edge_ids]
130        n2  = Normals[vol_ids,2*edge_ids+1]
131
132        # Transfer these quantities to the boundary array
133        Stage.boundary_values[ids]  = Stage.edge_values[vol_ids,edge_ids]
134        Elev.boundary_values[ids]   = Elev.edge_values[vol_ids,edge_ids]
135        Height.boundary_values[ids] = Height.edge_values[vol_ids,edge_ids]
136
137        # Rotate and negate Momemtum
138        q1 = Xmom.edge_values[vol_ids,edge_ids]
139        q2 = Ymom.edge_values[vol_ids,edge_ids]
140
141        r1 = -q1*n1 - q2*n2
142        r2 = -q1*n2 + q2*n1
143
144        Xmom.boundary_values[ids] = n1*r1 - n2*r2
145        Ymom.boundary_values[ids] = n2*r1 + n1*r2
146
147        # Rotate and negate Velocity
148        q1 = Xvel.edge_values[vol_ids,edge_ids]
149        q2 = Yvel.edge_values[vol_ids,edge_ids]
150
151        r1 = q1*n1 + q2*n2
152        r2 = q1*n2 - q2*n1
153
154        Xvel.boundary_values[ids] = n1*r1 - n2*r2
155        Yvel.boundary_values[ids] = n2*r1 + n1*r2
156
157
158
159class Transmissive_momentum_set_stage_boundary(Boundary):
160    """Returns same momentum conserved quantities as
161    those present in its neighbour volume.
162    Sets stage by specifying a function f of time which may either be a
163    vector function or a scalar function
164
165    Example:
166
167    def waveform(t):
168        return sea_level + normalized_amplitude/cosh(t-25)**2
169
170    Bts = Transmissive_momentum_set_stage_boundary(domain, waveform)
171
172    Underlying domain must be specified when boundary is instantiated
173    """
174
175    def __init__(self, domain=None, function=None):
176        Boundary.__init__(self)
177        """ Instantiate a Transmissive_momentum_set_stage_boundary. """
178
179
180        if domain is None:
181            msg = 'Domain must be specified for this type boundary'
182            raise Exception, msg
183
184        if function is None:
185            msg = 'Function must be specified for this type boundary'
186            raise Exception, msg
187
188        self.domain = domain
189
190        if isinstance(function, (int, float)):
191            tmp = function
192            function = lambda t: tmp
193           
194        self.function = function
195
196
197    def __repr__(self):
198        """ Return a representation of this instance. """
199        return 'Transmissive_momentum_set_stage_boundary(%s)' % self.domain
200
201    def evaluate(self, vol_id, edge_id):
202        """Transmissive momentum set stage boundaries return the edge momentum
203        values of the volume they serve.
204
205        vol_id is volume id
206        edge_id is the edge within the volume
207        """
208
209        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
210        t = self.domain.get_time()
211
212        if hasattr(self.function, 'time'):
213            # Roll boundary over if time exceeds
214            while t > self.function.time[-1]:
215                msg = 'WARNING: domain time %.2f has exceeded' % t
216                msg += 'time provided in '
217                msg += 'transmissive_momentum_set_stage_boundary object.\n'
218                msg += 'I will continue, reusing the object from t==0'
219                log.critical(msg)
220                t -= self.function.time[-1]
221
222        value = self.function(t)
223        try:
224            x = float(value)
225        except:
226            x = float(value[0])
227
228        q[0] = x
229           
230        return q
231
232        # FIXME: Consider this (taken from File_boundary) to allow
233        # spatial variation
234        # if vol_id is not None and edge_id is not None:
235        #     i = self.boundary_indices[ vol_id, edge_id ]
236        #     return self.F(t, point_id = i)
237        # else:
238        #     return self.F(t)
239
240
241class Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(Boundary):
242    """Returns the same normal momentum as that
243    present in neighbour volume edge. Zero out the tangential momentum.
244    Sets stage by specifying a function f of time which may either be a
245    vector function or a scalar function
246
247    Example:
248
249    def waveform(t):
250        return sea_level + normalized_amplitude/cosh(t-25)**2
251
252    Bts = Transmissive_n_momentum_zero_t_momentum_set_stage_boundary\
253                            (domain, waveform)
254
255    Underlying domain must be specified when boundary is instantiated
256    """
257
258    def __init__(self, domain=None, function=None):
259        """ Instantiate a
260            Transmissive_n_momentum_zero_t_momentum_set_stage_boundary.
261            domain is the domain containing the boundary
262            function is the function to apply
263        """
264
265        Boundary.__init__(self)
266
267        if domain is None:
268            msg = 'Domain must be specified for this type boundary'
269            raise Exception, msg
270
271        if function is None:
272            msg = 'Function must be specified for this type boundary'
273            raise Exception, msg
274
275        self.domain = domain
276        self.function = function
277
278
279    def __repr__(self):
280        """ Return a representation of this instance. """
281        msg = 'Transmissive_n_momentum_zero_t_momentum_set_stage_boundary'
282        msg += '(%s)' % self.domain
283        return msg
284
285
286    def evaluate(self, vol_id, edge_id):
287        """Transmissive_n_momentum_zero_t_momentum_set_stage_boundary
288        return the edge momentum values of the volume they serve.
289        """
290
291        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
292
293        normal = self.domain.get_normal(vol_id, edge_id)
294
295
296        t = self.domain.get_time()
297
298        if hasattr(self.function, 'time'):
299            # Roll boundary over if time exceeds
300            while t > self.function.time[-1]:
301                msg = 'WARNING: domain time %.2f has exceeded' % t
302                msg += 'time provided in '
303                msg += 'transmissive_momentum_set_stage_boundary object.\n'
304                msg += 'I will continue, reusing the object from t==0'
305                log.critical(msg)
306                t -= self.function.time[-1]
307
308        value = self.function(t)
309        try:
310            x = float(value)
311        except:
312            x = float(value[0])
313
314        q[0] = x
315        ndotq = (normal[0]*q[1] + normal[1]*q[2])
316        q[1] = normal[0]*ndotq
317        q[2] = normal[1]*ndotq
318
319        return q
320
321class Transmissive_stage_zero_momentum_boundary(Boundary):
322    """Return same stage as those present in its neighbour volume.
323    Set momentum to zero.
324
325    Underlying domain must be specified when boundary is instantiated
326    """
327
328    def __init__(self, domain=None):
329        """ Instantiate a Transmissive (zero momentum) boundary. """
330
331        Boundary.__init__(self)
332
333        if domain is None:
334            msg = ('Domain must be specified for '
335                   'Transmissive_stage_zero_momentum boundary')
336            raise Exception, msg
337
338        self.domain = domain
339
340    def __repr__(self):
341        """ Return a representation of this instance. """
342        return 'Transmissive_stage_zero_momentum_boundary(%s)' % self.domain
343
344
345    def evaluate(self, vol_id, edge_id):
346        """Calculate transmissive (zero momentum) results. """
347
348        q = self.domain.get_conserved_quantities(vol_id, edge=edge_id)
349
350        q[1] = q[2] = 0.0
351        return q
352
353
354
355class Time_stage_zero_momentum_boundary(Boundary):
356    """Time dependent boundary returns values for stage
357    conserved quantities as a function of time.
358    Must specify domain to get access to model time and a function of t
359    which must return conserved stage quantities as a function time.
360   
361    Example:
362      B = Time_stage_zero_momentum_boundary(domain,
363                        function=lambda t: (60<t<3660)*2)
364     
365      This will produce a boundary condition with is a 2m high square wave
366      starting 60 seconds into the simulation and lasting one hour.
367      Momentum applied will be 0 at all times.
368                       
369    """
370
371    def __init__(self, domain=None,
372                 #f=None, # Should be removed and replaced by function below
373                 function=None,
374                 default_boundary=None,
375                 verbose=False):
376        Boundary.__init__(self)
377
378        self.default_boundary = default_boundary
379        self.default_boundary_invoked = False    # Flag
380        self.domain = domain
381        self.verbose = verbose
382
383        if domain is None:
384            raise Exception('You must specify a domain to Time_stage_zero_momemtum_boundary')
385           
386        if function is None:
387            raise Exception('You must specify a function to Time_stage_zero_momemtum_boundary')
388           
389       
390        try:
391            q = function(0.0)
392        except Exception, e:
393            msg = 'Function for time stage boundary could not be executed:\n%s' %e
394            raise Exception(msg)
395
396
397        try:
398            q = float(q)
399        except:
400            msg = 'Return value from time boundary function could '
401            msg += 'not be converted into a float.\n'
402            msg += 'I got %s' %str(q)
403            raise Exception(msg)
404
405
406        self.f = function
407        self.domain = domain
408
409    def __repr__(self):
410        return 'Time_stage_zero_momemtum_boundary'
411
412    def get_time(self):
413
414        return self.domain.get_time()
415
416    def evaluate(self, vol_id=None, edge_id=None):
417
418        return self.get_boundary_values()
419
420
421    def evaluate_segment(self, domain, segment_edges):
422
423        if segment_edges is None:
424            return
425        if domain is None:
426            return
427
428        ids = segment_edges
429
430        vol_ids  = domain.boundary_cells[ids]
431        edge_ids = domain.boundary_edges[ids]
432
433        q_bdry = self.get_boundary_values()
434
435        #-------------------------------------------------
436        # Now update boundary values
437        #-------------------------------------------------
438        domain.quantities['stage'].boundary_values[ids] = q_bdry
439        domain.quantities['xmomentum'].boundary_values[ids] = 0.0
440        domain.quantities['ymomentum'].boundary_values[ids] = 0.0       
441
442
443    def get_boundary_values(self, t=None):
444
445        if t is None:
446            t = self.get_time()
447           
448        try:
449            res = self.f(t)
450        except Modeltime_too_early, e:
451            raise Modeltime_too_early(e)
452        except Modeltime_too_late, e:
453            if self.default_boundary is None:
454                raise Exception(e) # Reraise exception
455            else:
456                # Pass control to default boundary
457                res = self.default_boundary
458               
459                # Ensure that result cannot be manipulated
460                # This is a real danger in case the
461                # default_boundary is a Dirichlet type
462                # for instance.
463                res = res.copy() 
464               
465                if self.default_boundary_invoked is False:
466                    if self.verbose:               
467                        # Issue warning the first time
468                        msg = '%s' %str(e)
469                        msg += 'Instead I will use the default boundary: %s\n'\
470                            %str(self.default_boundary) 
471                        msg += 'Note: Further warnings will be suppressed'
472                        log.critical(msg)
473               
474                    # FIXME (Ole): Replace this crude flag with
475                    # Python's ability to print warnings only once.
476                    # See http://docs.python.org/lib/warning-filter.html
477                    self.default_boundary_invoked = True
478
479        return res
480
481
482class Characteristic_stage_boundary(Boundary):
483    """Sets the stage via a function and the momentum is determined
484    via assumption of simple incoming wave (uses Riemann invariant)
485
486
487    Example:
488
489    def waveform(t):
490        return sea_level + normalized_amplitude/cosh(t-25)**2
491
492    Bcs = Characteristic_stage_boundary(domain, waveform)
493
494    Underlying domain must be specified when boundary is instantiated
495    """
496
497    def __init__(self, domain=None, function=None, default_stage = 0.0):
498        """ Instantiate a
499            Characteristic_stage_boundary.
500            domain is the domain containing the boundary
501            function is the function to apply the wave
502            default_stage is the assumed stage pre the application of wave
503        """
504
505        Boundary.__init__(self)
506
507        if domain is None:
508            msg = 'Domain must be specified for this type boundary'
509            raise Exception, msg
510
511        if function is None:
512            msg = 'Function must be specified for this type boundary'
513            raise Exception, msg
514
515        self.domain = domain
516        self.function = function
517        self.default_stage = default_stage
518
519        self.Elev  = domain.quantities['elevation']
520        self.Stage = domain.quantities['stage']
521        self.Height = domain.quantities['height']
522
523    def __repr__(self):
524        """ Return a representation of this instance. """
525        msg = 'Characteristic_stage_boundary '
526        msg += '(%s) ' % self.domain
527        msg += '(%s) ' % self.default_stage
528        return msg
529
530
531    def evaluate(self, vol_id, edge_id):
532        """Calculate reflections (reverse outward momentum).
533
534        vol_id   
535        edge_id 
536        """
537       
538        t = self.domain.get_time()
539
540
541        value = self.function(t)
542        try:
543            stage = float(value)
544        except:
545            stage = float(value[0])
546
547
548
549        q = self.conserved_quantities
550        #q[0] = self.stage[vol_id, edge_id]
551        q[0] = stage
552        q[1] = self.xmom[vol_id, edge_id]
553        q[2] = self.ymom[vol_id, edge_id]
554
555        normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
556
557        r = rotate(q, normal, direction = 1)
558        r[1] = -r[1]
559        q = rotate(r, normal, direction = -1)
560
561
562        return q
563
564
565
566
567
568
569    def evaluate_segment(self, domain, segment_edges):
570        """Apply reflective BC on the boundary edges defined by
571        segment_edges
572        """
573
574        if segment_edges is None:
575            return
576        if domain is None:
577            return
578
579        t = self.domain.get_time()
580
581
582        value = self.function(t)
583        try:
584            stage = float(value)
585        except:
586            stage = float(value[0])
587                       
588
589        ids = segment_edges
590        vol_ids  = domain.boundary_cells[ids]
591        edge_ids = domain.boundary_edges[ids]
592
593        Stage = domain.quantities['stage']
594        Elev  = domain.quantities['elevation']
595        Height= domain.quantities['height']
596        Xmom  = domain.quantities['xmomentum']
597        Ymom  = domain.quantities['ymomentum']
598        Xvel  = domain.quantities['xvelocity']
599        Yvel  = domain.quantities['yvelocity']
600
601        Normals = domain.normals
602
603        #print vol_ids
604        #print edge_ids
605        #Normals.reshape((4,3,2))
606        #print Normals.shape
607        #print Normals[vol_ids, 2*edge_ids]
608        #print Normals[vol_ids, 2*edge_ids+1]
609       
610        n1  = Normals[vol_ids,2*edge_ids]
611        n2  = Normals[vol_ids,2*edge_ids+1]
612
613        # Transfer these quantities to the boundary array
614        Stage.boundary_values[ids]  = Stage.edge_values[vol_ids,edge_ids]
615        Elev.boundary_values[ids]   = Elev.edge_values[vol_ids,edge_ids]
616        Height.boundary_values[ids] = Height.edge_values[vol_ids,edge_ids]
617
618        # Rotate and negate Momemtum
619        q1 = Xmom.edge_values[vol_ids,edge_ids]
620        q2 = Ymom.edge_values[vol_ids,edge_ids]
621
622        r1 = -q1*n1 - q2*n2
623        r2 = -q1*n2 + q2*n1
624
625        Xmom.boundary_values[ids] = n1*r1 - n2*r2
626        Ymom.boundary_values[ids] = n2*r1 + n1*r2
627
628        # Rotate and negate Velocity
629        q1 = Xvel.edge_values[vol_ids,edge_ids]
630        q2 = Yvel.edge_values[vol_ids,edge_ids]
631
632        r1 = q1*n1 + q2*n2
633        r2 = q1*n2 - q2*n1
634
635        Xvel.boundary_values[ids] = n1*r1 - n2*r2
636        Yvel.boundary_values[ids] = n2*r1 + n1*r2
637
638       
639
640
641class Dirichlet_discharge_boundary(Boundary):
642    """ Class for a Dirichlet discharge boundary.
643
644    Sets stage (stage0)
645    Sets momentum (wh0) in the inward normal direction.
646
647    Underlying domain must be specified when boundary is instantiated
648    """
649
650    def __init__(self, domain=None, stage0=None, wh0=None):
651        Boundary.__init__(self)
652        """ Instantiate a Dirichlet discharge boundary.
653            domain underlying domain
654            stage0 stag
655            wh0 momentum in the inward normal direction.
656            """
657
658        if domain is None:
659            msg = 'Domain must be specified for this type of boundary'
660            raise Exception, msg
661
662        if stage0 is None:
663            raise Exception, 'Stage must be specified for this type of boundary'
664
665        if wh0 is None:
666            wh0 = 0.0
667
668        self.domain = domain
669        self.stage0 = stage0
670        self.wh0 = wh0
671
672
673    def __repr__(self):
674        """ Return a representation of this instance. """
675        return 'Dirichlet_Discharge_boundary(%s)' % self.domain
676
677    def evaluate(self, vol_id, edge_id):
678        """Set discharge in the (inward) normal direction"""
679
680        normal = self.domain.get_normal(vol_id,edge_id)
681        q = [self.stage0, -self.wh0*normal[0], -self.wh0*normal[1]]
682        return q
683
684        # FIXME: Consider this (taken from File_boundary) to allow
685        # spatial variation
686        # if vol_id is not None and edge_id is not None:
687        #     i = self.boundary_indices[ vol_id, edge_id ]
688        #     return self.F(t, point_id = i)
689        # else:
690        #     return self.F(t)
691
692
693class Inflow_boundary(Boundary):
694    """Apply given flow in m^3/s to boundary segment.
695    Depth and momentum is derived using Manning's formula.
696
697    Underlying domain must be specified when boundary is instantiated
698    """
699   
700    # FIXME (Ole): This is work in progress and definitely not finished.
701    # The associated test has been disabled
702
703    def __init__(self, domain=None, rate=0.0):
704        Boundary.__init__(self)
705
706        if domain is None:
707            msg = 'Domain must be specified for '
708            msg += 'Inflow boundary'
709            raise Exception, msg
710
711        self.domain = domain
712       
713        # FIXME(Ole): Allow rate to be time dependent as well
714        self.rate = rate
715        self.tag = None # Placeholder for tag associated with this object.
716
717    def __repr__(self):
718        return 'Inflow_boundary(%s)' %self.domain
719
720    def evaluate(self, vol_id, edge_id):
721        """Apply inflow rate at each edge of this boundary
722        """
723       
724        # First find all segments having the same tag is vol_id, edge_id
725        # This will be done the first time evaluate is called.
726        if self.tag is None:
727            boundary = self.domain.boundary
728            self.tag = boundary[(vol_id, edge_id)]       
729           
730            # Find total length of boundary with this tag
731            length = 0.0
732            for v_id, e_id in boundary:
733                if self.tag == boundary[(v_id, e_id)]:
734                    length += self.domain.mesh.get_edgelength(v_id, e_id)           
735
736            self.length = length
737            self.average_momentum = self.rate/length
738           
739           
740        # Average momentum has now been established across this boundary
741        # Compute momentum in the inward normal direction
742       
743        inward_normal = -self.domain.mesh.get_normal(vol_id, edge_id)       
744        xmomentum, ymomentum = self.average_momentum * inward_normal
745           
746        # Compute depth based on Manning's formula v = 1/n h^{2/3} sqrt(S)
747        # Where v is velocity, n is manning's coefficient, h is depth
748        # and S is the slope into the domain.
749        # Let mu be the momentum (vh), then this equation becomes:
750        #            mu = 1/n h^{5/3} sqrt(S)
751        # from which we can isolate depth to get
752        #             h = (mu n/sqrt(S) )^{3/5}
753       
754        slope = 0 # get gradient for this triangle dot normal
755       
756        # get manning coef from this triangle
757        friction = self.domain.get_quantity('friction').get_values(\
758                    location='edges', indices=[vol_id])[0]
759        mannings_n = friction[edge_id]
760
761        if slope > epsilon and mannings_n > epsilon:
762            depth = pow(self.average_momentum * mannings_n/math.sqrt(slope), \
763                        3.0/5) 
764        else:
765            depth = 1.0
766           
767        # Elevation on this edge   
768       
769        z = self.domain.get_quantity('elevation').get_values(\
770                    location='edges', indices=[vol_id])[0]
771        elevation = z[edge_id]
772           
773        # Assign conserved quantities and return
774        q = num.array([elevation + depth, xmomentum, ymomentum], num.float)
775        return q
776
777
778       
779   
780           
781       
782class Field_boundary(Boundary):
783    """Set boundary from given field represented in an sww file containing
784    values for stage, xmomentum and ymomentum.
785
786    Optionally, the user can specify mean_stage to offset the stage provided
787    in the sww file.
788
789    This function is a thin wrapper around the generic File_boundary. The
790    difference between the File_boundary and Field_boundary is only that the
791    Field_boundary will allow you to change the level of the stage height when
792    you read in the boundary condition. This is very useful when running
793    different tide heights in the same area as you need only to convert one
794    boundary condition to a SWW file, ideally for tide height of 0 m
795    (saving disk space). Then you can use Field_boundary to read this SWW file
796    and change the stage height (tide) on the fly depending on the scenario.
797    """
798
799    def __init__(self,
800                 filename,
801                 domain,
802                 mean_stage=0.0,
803                 time_thinning=1,
804                 time_limit=None,
805                 boundary_polygon=None,
806                 default_boundary=None,
807                 use_cache=False,
808                 verbose=False):
809        """Constructor
810
811        filename: Name of sww file containing stage and x/ymomentum
812        domain: pointer to shallow water domain for which the boundary applies
813        mean_stage: The mean water level which will be added to stage derived
814                    from the boundary condition
815        time_thinning: Will set how many time steps from the sww file read in
816                       will be interpolated to the boundary. For example if
817                       the sww file has 1 second time steps and is 24 hours
818                       in length it has 86400 time steps. If you set
819                       time_thinning to 1 it will read all these steps.
820                       If you set it to 100 it will read every 100th step eg
821                       only 864 step. This parameter is very useful to increase
822                       the speed of a model run that you are setting up
823                       and testing.
824
825        default_boundary: Must be either None or an instance of a
826                          class descending from class Boundary.
827                          This will be used in case model time exceeds
828                          that available in the underlying data.
829
830                          Note that mean_stage will also be added to this.
831
832        time_limit:
833        boundary_polygon:
834        use_cache:        True if caching is to be used.
835        verbose:          True if this method is to be verbose.
836
837        """
838
839        # Create generic file_boundary object
840        self.file_boundary = File_boundary(filename,
841                                           domain,
842                                           time_thinning=time_thinning,
843                                           time_limit=time_limit,
844                                           boundary_polygon=boundary_polygon,
845                                           default_boundary=default_boundary,
846                                           use_cache=use_cache,
847                                           verbose=verbose)
848
849        # Record information from File_boundary
850        self.F = self.file_boundary.F
851        self.domain = self.file_boundary.domain
852
853        # Record mean stage
854        self.mean_stage = mean_stage
855
856
857    def __repr__(self):
858        """ Generate a string representation of this instance. """
859        return 'Field boundary'
860
861
862    def evaluate(self, vol_id=None, edge_id=None):
863        """ Calculate 'field' boundary results.
864            vol_id and edge_id are ignored
865
866            Return linearly interpolated values based on domain.time
867        """
868
869        # Evaluate file boundary
870        q = self.file_boundary.evaluate(vol_id, edge_id)
871
872        # Adjust stage
873        for j, name in enumerate(self.domain.conserved_quantities):
874            if name == 'stage':
875                q[j] += self.mean_stage
876        return q
877
878
879
880
881
882class Flather_external_stage_zero_velocity_boundary(Boundary):
883    """
884
885    Boundary condition based on a Flather type approach
886    (setting the external stage with a function, and a zero external velocity),
887   
888
889    The idea is similar (but not identical) to that described on page 239 of
890    the following article:
891
892    Article{blayo05,
893      Title                    = {Revisiting open boundary conditions from the point of view of characteristic variables},
894      Author                   = {Blayo, E. and Debreu, L.},
895      Journal                  = {Ocean Modelling},
896      Year                     = {2005},
897      Pages                    = {231-252},
898      Volume                   = {9},
899    }
900
901    Approach
902    1) The external (outside boundary) stage is set with a function, the
903       external velocity is zero, the internal stage and velocity are taken from the
904       domain values.
905    2) Some 'characteristic like' variables are computed, depending on whether
906       the flow is incoming or outgoing. See Blayo and Debreu (2005)
907    3) The boundary conserved quantities are computed from these characteristic
908       like variables
909
910    This has been useful as a 'weakly reflecting' boundary when the stage should
911    be approximately specified but allowed to adapt to outgoing waves.
912
913
914    Example:
915
916    def waveform(t):
917        return sea_level + normalized_amplitude/cosh(t-25)**2
918
919    Bf = Flather_external_stage_zero_velocity_boundary(domain, waveform)
920
921    Underlying domain must be specified when boundary is instantiated
922
923
924   
925    """
926
927    def __init__(self, domain=None, function=None):
928        """ Instantiate a
929            Nudge_boundary.
930            domain is the domain containing the boundary
931            function is the function to apply
932        """
933
934        Boundary.__init__(self)
935
936        if domain is None:
937            msg = 'Domain must be specified for this type boundary'
938            raise Exception, msg
939
940        if function is None:
941            msg = 'Function must be specified for this type boundary'
942            raise Exception, msg
943
944        self.domain = domain
945        self.function = function
946
947
948    def __repr__(self):
949        """ Return a representation of this instance. """
950        msg = 'Nudge_boundary'
951        msg += '(%s)' % self.domain
952        return msg
953
954
955    def evaluate(self, vol_id, edge_id):
956        """
957        """
958
959        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
960        bed = self.domain.quantities['elevation'].centroid_values[vol_id]
961        depth_inside=max(q[0]-bed,0.0)
962        dt=self.domain.timestep
963
964        normal = self.domain.get_normal(vol_id, edge_id)
965
966
967        t = self.domain.get_time()
968
969        value = self.function(t)
970        try:
971            stage_outside = float(value)
972        except:
973            stage_outside = float(value[0])
974
975        if(depth_inside==0.):
976            q[0] = stage_outside
977            q[1] = 0. 
978            q[2] = 0. 
979
980        else:
981
982            # Asssume sub-critical flow. Set the values of the characteristics as
983            # appropriate, depending on whether we have inflow or outflow
984
985            # These calculations are based on the paper cited above
986            sqrt_g_on_depth_inside = (gravity/depth_inside)**0.5
987            ndotq_inside = (normal[0]*q[1] + normal[1]*q[2]) # momentum perpendicular to the boundary
988            if(ndotq_inside>0.):
989                # Outflow (assumed subcritical)
990                # Compute characteristics using a particular extrapolation
991                #
992                # Theory: 2 characteristics coming from inside domain, only
993                # need to impose one characteristic from outside
994                #
995
996                # w1 =  u - sqrt(g/depth)*(Stage_outside)  -- uses 'outside' info
997                w1 = 0. - sqrt_g_on_depth_inside*stage_outside
998
999                # w2 = v [velocity parallel to boundary] -- uses 'inside' info
1000                w2 = (+normal[1]*q[1] -normal[0]*q[2])/depth_inside
1001
1002                # w3 = u + sqrt(g/depth)*(Stage_inside) -- uses 'inside info'
1003                w3 = ndotq_inside/depth_inside + sqrt_g_on_depth_inside*q[0]
1004               
1005            else:
1006                # Inflow (assumed subcritical)
1007                # Need to set 2 characteristics from outside information
1008               
1009                # w1 =  u - sqrt(g/depth)*(Stage_outside)  -- uses 'outside' info
1010                w1 = 0. - sqrt_g_on_depth_inside*stage_outside
1011
1012                # w2 = v [velocity parallel to boundary] -- uses 'outside' info
1013                w2 = 0.
1014
1015                # w3 = u + sqrt(g/depth)*(Stage_inside) -- uses 'inside info'
1016                w3 = ndotq_inside/depth_inside + sqrt_g_on_depth_inside*q[0]
1017
1018
1019            q[0] = (w3-w1)/(2*sqrt_g_on_depth_inside)
1020            qperp= (w3+w1)/2.*depth_inside
1021            qpar=  w2*depth_inside
1022
1023            # So q[1], q[2] = qperp*(normal[0], normal[1]) + qpar*(-normal[1], normal[0])
1024
1025            q[1] = qperp*normal[0] + qpar*normal[1]
1026            q[2] = qperp*normal[1] -qpar*normal[0]
1027
1028        return q
Note: See TracBrowser for help on using the repository browser.