source: anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py @ 7573

Last change on this file since 7573 was 7573, checked in by steve, 14 years ago

Committing a version of shallow_water_balanced which passes it unit test
using a version of edge limiting which doesn't limit boundary edges. THis
is useful to allow linear functions to be reconstructed.

Had to play with the transmissive BC and use centroid values which is
set via domain set_centroid_transmissive_bc

File size: 24.7 KB
Line 
1"""boundary.py - Classes for implementing boundary conditions
2"""
3
4from warnings import warn
5
6from anuga.utilities.numerical_tools import NAN   
7from anuga.fit_interpolate.interpolate import Modeltime_too_late
8from anuga.fit_interpolate.interpolate import Modeltime_too_early
9import anuga.utilities.log as log
10
11import numpy as num
12
13
14class Boundary:
15    """Base class for boundary conditions.
16       Specific boundary conditions must provide values for
17       the conserved_quantities
18
19       A boundary object has one neighbour; the one it
20       serves.
21    """
22
23    def __init__(self):
24        pass
25
26    def evaluate(self, vol_id=None, edge_id=None):
27        msg = 'Generic class Boundary must be subclassed'
28        raise Exception, msg
29
30
31class Transmissive_boundary(Boundary):
32    """Transmissive boundary returns same conserved quantities as
33    those present in its neighbour volume.
34
35    Underlying domain must be specified when boundary is instantiated
36    """
37
38    def __init__(self, domain = None):
39        Boundary.__init__(self)
40
41        if domain is None:
42            msg = 'Domain must be specified for transmissive boundary'
43            raise Exception, msg
44
45        self.domain = domain
46
47    def __repr__(self):
48        return 'Transmissive_boundary(%s)' %self.domain
49
50    def evaluate(self, vol_id, edge_id):
51        """Transmissive boundaries return the edge values
52        of the volume they serve.
53        """
54
55        if self.domain.centroid_transmissive_bc :
56            q = self.domain.get_conserved_quantities(vol_id)
57        else:
58            q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
59
60
61        return q
62
63
64class Dirichlet_boundary(Boundary):
65    """Dirichlet boundary returns constant values for the
66    conserved quantities
67    """
68
69
70    def __init__(self, conserved_quantities=None):
71        Boundary.__init__(self)
72
73        if conserved_quantities is None:
74            msg = 'Must specify one value for each conserved quantity'
75            raise Exception, msg
76
77        self.conserved_quantities=num.array(conserved_quantities, num.float)
78
79    def __repr__(self):
80        return 'Dirichlet boundary (%s)' %self.conserved_quantities
81
82    def evaluate(self, vol_id=None, edge_id=None):
83        return self.conserved_quantities
84
85
86
87class Time_boundary(Boundary):
88    """Time dependent boundary returns values for the
89    conserved quantities as a function of time.
90    Must specify domain to get access to model time and a function of t
91    which must return conserved quantities as a function time.
92   
93    Example:
94      B = Time_boundary(domain,
95                        function=lambda t: [(60<t<3660)*2, 0, 0])
96     
97      This will produce a boundary condition with is a 2m high square wave
98      starting 60 seconds into the simulation and lasting one hour.
99      Momentum applied will be 0 at all times.
100                       
101    """
102
103    # FIXME (Ole): We should rename f to function to be consistent with
104    # Transmissive_Momentum_Set_Stage_Boundary (cf posting by rrraman)
105    def __init__(self, domain=None,
106                 f=None, # Should be removed and replaced by function below
107                 function=None,
108                 default_boundary=None,
109                 verbose=False):
110        Boundary.__init__(self)
111
112        self.default_boundary = default_boundary
113        self.default_boundary_invoked = False    # Flag
114        self.domain = domain
115        self.verbose = verbose
116
117        if domain is None:
118            raise Exception, 'You must specify a domain to Time_boundary'
119
120        # FIXME: Temporary code to deal with both f and function
121        if function is not None and f is not None:
122            raise Exception, 'Specify either function or f to Time_boundary'
123           
124        if function is None and f is None:
125            raise Exception, 'You must specify a function to Time_boundary'
126           
127        if f is None:
128            f = function
129        #####
130       
131        try:
132            q = f(0.0)
133        except Exception, e:
134            msg = 'Function for time boundary could not be executed:\n%s' %e
135            raise msg
136
137
138        try:
139            q = num.array(q, num.float)
140        except:
141            msg = 'Return value from time boundary function could '
142            msg += 'not be converted into a numeric array of floats.\n'
143            msg += 'Specified function should return either list or array.\n'
144            msg += 'I got %s' %str(q)
145            raise msg
146
147        msg = 'ERROR: Time boundary function must return a 1d list or array '
148        assert len(q.shape) == 1, msg
149
150        d = len(domain.conserved_quantities)
151        msg = 'Return value for function must be a list or an array of length %d' %d
152        assert len(q) == d, msg
153
154        self.f = f
155        self.domain = domain
156
157    def __repr__(self):
158        return 'Time boundary'
159
160    def evaluate(self, vol_id=None, edge_id=None):
161        # FIXME (Ole): I think this should be get_time(), see ticket:306
162        try:
163            res = self.f(self.domain.time)
164        except Modeltime_too_early, e:
165            raise Modeltime_too_early, e
166        except Modeltime_too_late, e:
167            if self.default_boundary is None:
168                raise Exception, e # Reraise exception
169            else:
170                # Pass control to default boundary
171                res = self.default_boundary.evaluate(vol_id, edge_id)
172               
173                # Ensure that result cannot be manipulated
174                # This is a real danger in case the
175                # default_boundary is a Dirichlet type
176                # for instance.
177                res = res.copy() 
178               
179                if self.default_boundary_invoked is False:
180                    if self.verbose:               
181                        # Issue warning the first time
182                        msg = '%s' %str(e)
183                        msg += 'Instead I will use the default boundary: %s\n'\
184                            %str(self.default_boundary) 
185                        msg += 'Note: Further warnings will be supressed'
186                        log.critical(msg)
187               
188                    # FIXME (Ole): Replace this crude flag with
189                    # Python's ability to print warnings only once.
190                    # See http://docs.python.org/lib/warning-filter.html
191                    self.default_boundary_invoked = True
192
193        return res
194
195class Time_space_boundary(Boundary):
196    """Time and spatially dependent boundary returns values for the
197    conserved quantities as a function of time and space.
198    Must specify domain to get access to model time and a function of t,x,y
199    which must return conserved quantities at specified time and location.
200   
201    Example:
202      B = Time_space_boundary(domain,
203                              function=lambda t,x,y: [(60<t<3660)*2, 0, 0])
204     
205      This will produce a boundary condition with is a 2m high square wave
206      starting 60 seconds into the simulation and lasting one hour.
207      Momentum applied will be 0 at all times.
208                       
209    """
210
211
212    def __init__(self, domain=None,
213                 function=None,
214                 default_boundary=None,
215                 verbose=False):
216        Boundary.__init__(self)
217
218        self.default_boundary = default_boundary
219        self.default_boundary_invoked = False    # Flag
220        self.domain = domain
221        self.verbose = verbose
222
223       
224        if function is None:
225            raise Exception, 'You must specify a function to Time_space_boundary'
226           
227        try:
228            q = function(0.0, 0.0, 0.0)
229        except Exception, e:
230            msg = 'Function for time_space_boundary could not be executed:\n%s' %e
231            raise msg
232
233
234        try:
235            q = num.array(q, num.float)
236        except:
237            msg = 'Return value from time_space_boundary function could '
238            msg += 'not be converted into a numeric array of floats.\n'
239            msg += 'Specified function should return either list or array.\n'
240            msg += 'I got %s' %str(q)
241            raise msg
242
243        msg = 'ERROR: Time_space_boundary function must return a 1d list or array '
244        assert len(q.shape) == 1, msg
245
246        d = len(domain.conserved_quantities)
247        msg = 'Return value for function must be a list or an array of length %d' %d
248        assert len(q) == d, msg
249
250        self.function = function
251        self.domain = domain
252
253    def __repr__(self):
254        return 'Time space boundary'
255
256    def evaluate(self, vol_id=None, edge_id=None):
257        # FIXME (Ole): I think this should be get_time(), see ticket:306
258
259        [x,y] = self.domain.get_edge_midpoint_coordinate(vol_id, edge_id)
260
261        try:
262            res = self.function(self.domain.get_time(), x, y)
263        except Modeltime_too_early, e:
264            raise Modeltime_too_early, e
265        except Modeltime_too_late, e:
266            if self.default_boundary is None:
267                raise Exception, e # Reraise exception
268            else:
269                # Pass control to default boundary
270                res = self.default_boundary.evaluate(vol_id, edge_id)
271               
272                # Ensure that result cannot be manipulated
273                # This is a real danger in case the
274                # default_boundary is a Dirichlet type
275                # for instance.
276                res = res.copy() 
277               
278                if self.default_boundary_invoked is False:
279                    if self.verbose:               
280                        # Issue warning the first time
281                        msg = '%s' %str(e)
282                        msg += 'Instead I will use the default boundary: %s\n'\
283                            %str(self.default_boundary) 
284                        msg += 'Note: Further warnings will be supressed'
285                        log.critical(msg)
286               
287                    # FIXME (Ole): Replace this crude flag with
288                    # Python's ability to print warnings only once.
289                    # See http://docs.python.org/lib/warning-filter.html
290                    self.default_boundary_invoked = True
291
292        return res
293
294class File_boundary(Boundary):
295    """The File_boundary reads values for the conserved
296    quantities from an sww NetCDF file, and returns interpolated values
297    at the midpoints of each associated boundary segment.
298    Time dependency is interpolated linearly.
299
300    Assumes that file contains a time series and possibly
301    also spatial info. See docstring for File_function in util.py
302    for details about admissible file formats
303
304    File boundary must read and interpolate from *smoothed* version
305    as stored in sww and cannot work with the discontinuos triangles.
306
307    Example:
308    Bf = File_boundary('source_file.sww', domain)
309
310
311    Note that the resulting solution history is not exactly the same as if
312    the models were coupled as there is no feedback into the source model.
313   
314    Optional keyword argument default_boundary must be either None or
315    an instance of class descending from class Boundary.
316    This will be used in case model time exceeds that available in the
317    underlying data.
318       
319    """
320
321    # FIXME (Ole): I kind of like the name Spatio_Temporal_boundary for this
322    # rather than File_boundary
323
324    def __init__(self, filename, domain, 
325                 time_thinning=1, 
326                 time_limit=None,
327                 boundary_polygon=None,   
328                 default_boundary=None,
329                 use_cache=False, 
330                 verbose=False): 
331
332        import time
333        from anuga.config import time_format
334        from anuga.abstract_2d_finite_volumes.util import file_function
335
336        Boundary.__init__(self)
337
338        # Get x,y vertex coordinates for all triangles
339        V = domain.vertex_coordinates
340
341        # Compute midpoint coordinates for all boundary elements
342        # Only a subset may be invoked when boundary is evaluated but
343        # we don't know which ones at this stage since this object can
344        # be attached to
345        # any tagged boundary later on.
346
347        if verbose: log.critical('Find midpoint coordinates of entire boundary')
348        self.midpoint_coordinates = num.zeros((len(domain.boundary), 2), num.float)
349        boundary_keys = domain.boundary.keys()
350
351        xllcorner = domain.geo_reference.get_xllcorner()
352        yllcorner = domain.geo_reference.get_yllcorner()       
353       
354
355        # Make ordering unique #FIXME: should this happen in domain.py?
356        boundary_keys.sort()
357
358        # Record ordering #FIXME: should this also happen in domain.py or general_mesh.py?
359        self.boundary_indices = {}
360        for i, (vol_id, edge_id) in enumerate(boundary_keys):
361
362            self.midpoint_coordinates[i,:] = domain.get_edge_midpoint_coordinate(vol_id, edge_id, 
363                                                                                 absolute=True)
364
365            # Register index of this boundary edge for use with evaluate
366            self.boundary_indices[(vol_id, edge_id)] = i
367
368           
369           
370        if verbose: log.critical('Initialise file_function')
371        self.F = file_function(filename,
372                               domain,
373                               quantities=domain.conserved_quantities,
374                               interpolation_points=self.midpoint_coordinates,
375                               time_thinning=time_thinning,
376                               time_limit=time_limit,
377                               use_cache=use_cache, 
378                               verbose=verbose,
379                               boundary_polygon=boundary_polygon)
380                             
381        # Check and store default_boundary
382        msg = 'Keyword argument default_boundary must be either None '
383        msg += 'or a boundary object.\n I got %s' %(str(default_boundary))
384        assert default_boundary is None or\
385            isinstance(default_boundary, Boundary), msg
386        self.default_boundary = default_boundary
387        self.default_boundary_invoked = False    # Flag
388
389        # Store pointer to domain and verbosity
390        self.domain = domain
391        self.verbose = verbose
392
393
394        # Here we'll flag indices outside the mesh as a warning
395        # as suggested by Joaquim Luis in sourceforge posting
396        # November 2007
397        # We won't make it an error as it is conceivable that
398        # only part of mesh boundary is actually used with a given
399        # file boundary sww file.
400        if hasattr(self.F, 'indices_outside_mesh') and\
401               len(self.F.indices_outside_mesh) > 0:
402            msg = 'WARNING: File_boundary has points outside the mesh '
403            msg += 'given in %s. ' %filename
404            msg += 'See warning message issued by Interpolation_function '
405            msg += 'for details (should appear above somewhere if '
406            msg += 'verbose is True).\n'
407            msg += 'This is perfectly OK as long as the points that are '
408            msg += 'outside aren\'t used on the actual boundary segment.'
409            if verbose is True:           
410                log.critical(msg)
411            #raise Exception(msg)
412
413        # Test that file function can be called
414        q = self.F(0, point_id=0)
415        d = len(domain.conserved_quantities)
416        msg = 'Values specified in file %s must be ' %filename
417        msg += ' a list or an array of length %d' %d
418        assert len(q) == d, msg
419
420
421    def __repr__(self):
422        return 'File boundary'
423
424
425    def evaluate(self, vol_id=None, edge_id=None):
426        """Return linearly interpolated values based on domain.time
427        at midpoint of segment defined by vol_id and edge_id.
428        """
429
430        # FIXME (Ole): I think this should be get_time(), see ticket:306
431        t = self.domain.time
432       
433        if vol_id is not None and edge_id is not None:
434            i = self.boundary_indices[ vol_id, edge_id ]
435           
436            try:
437                res = self.F(t, point_id=i)
438            except Modeltime_too_early, e:
439                raise Modeltime_too_early, e
440            except Modeltime_too_late, e:
441                if self.default_boundary is None:
442                    raise Exception, e # Reraise exception
443                else:
444                    # Pass control to default boundary
445                    res = self.default_boundary.evaluate(vol_id, edge_id)
446                   
447                    # Ensure that result cannot be manipulated
448                    # This is a real danger in case the
449                    # default_boundary is a Dirichlet type
450                    # for instance.
451                    res = res.copy() 
452                   
453                    if self.default_boundary_invoked is False:
454                        # Issue warning the first time
455                        if self.verbose:
456                            msg = '%s' %str(e)
457                            msg += 'Instead I will use the default boundary: %s\n'\
458                                %str(self.default_boundary) 
459                            msg += 'Note: Further warnings will be supressed'
460                            log.critical(msg)
461                   
462                        # FIXME (Ole): Replace this crude flag with
463                        # Python's ability to print warnings only once.
464                        # See http://docs.python.org/lib/warning-filter.html
465                        self.default_boundary_invoked = True
466           
467            if num.any(res == NAN):
468                x,y=self.midpoint_coordinates[i,:]
469                msg = 'NAN value found in file_boundary at '
470                msg += 'point id #%d: (%.2f, %.2f).\n' %(i, x, y)
471
472                if hasattr(self.F, 'indices_outside_mesh') and\
473                       len(self.F.indices_outside_mesh) > 0:
474                    # Check if NAN point is due it being outside
475                    # boundary defined in sww file.
476
477                    if i in self.F.indices_outside_mesh:
478                        msg += 'This point refers to one outside the '
479                        msg += 'mesh defined by the file %s.\n'\
480                               %self.F.filename
481                        msg += 'Make sure that the file covers '
482                        msg += 'the boundary segment it is assigned to '
483                        msg += 'in set_boundary.'
484                    else:
485                        msg += 'This point is inside the mesh defined '
486                        msg += 'the file %s.\n' %self.F.filename
487                        msg += 'Check this file for NANs.'
488                raise Exception, msg
489           
490            return res
491        else:
492            msg = 'Boundary call without point_id not implemented.\n'
493            msg += 'vol_id=%s, edge_id=%s' %(str(vol_id), str(edge_id))
494            raise Exception, msg
495
496class AWI_boundary(Boundary):
497    """The AWI_boundary reads values for the conserved
498    quantities (only STAGE) from an sww NetCDF file, and returns interpolated values
499    at the midpoints of each associated boundary segment.
500    Time dependency is interpolated linearly.
501
502    Assumes that file contains a time series and possibly
503    also spatial info. See docstring for File_function in util.py
504    for details about admissible file formats
505
506    AWI boundary must read and interpolate from *smoothed* version
507    as stored in sww and cannot work with the discontinuos triangles.
508
509    Example:
510    Ba = AWI_boundary('source_file.sww', domain)
511
512
513    Note that the resulting solution history is not exactly the same as if
514    the models were coupled as there is no feedback into the source model.
515
516    This was added by Nils Goseberg et al in April 2009
517    """
518
519    def __init__(self, filename, domain, time_thinning=1, 
520                 use_cache=False, verbose=False):
521        import time
522        from anuga.config import time_format
523        from anuga.abstract_2d_finite_volumes.util import file_function
524
525        Boundary.__init__(self)
526
527        # Get x,y vertex coordinates for all triangles
528        V = domain.vertex_coordinates
529
530        # Compute midpoint coordinates for all boundary elements
531        # Only a subset may be invoked when boundary is evaluated but
532        # we don't know which ones at this stage since this object can
533        # be attached to
534        # any tagged boundary later on.
535
536        if verbose: log.critical('Find midpoint coordinates of entire boundary')
537        self.midpoint_coordinates = num.zeros((len(domain.boundary), 2), num.float)
538        boundary_keys = domain.boundary.keys()
539
540        xllcorner = domain.geo_reference.get_xllcorner()
541        yllcorner = domain.geo_reference.get_yllcorner()       
542
543        # Make ordering unique #FIXME: should this happen in domain.py?
544        boundary_keys.sort()
545
546        # Record ordering #FIXME: should this also happen in domain.py?
547        self.boundary_indices = {}
548        for i, (vol_id, edge_id) in enumerate(boundary_keys):
549            base_index = 3*vol_id
550            x0, y0 = V[base_index, :]
551            x1, y1 = V[base_index+1, :]
552            x2, y2 = V[base_index+2, :]
553           
554            # Compute midpoints
555            if edge_id == 0: m = num.array([(x1 + x2)/2, (y1 + y2)/2])
556            if edge_id == 1: m = num.array([(x0 + x2)/2, (y0 + y2)/2])
557            if edge_id == 2: m = num.array([(x1 + x0)/2, (y1 + y0)/2])
558
559            # Convert to absolute UTM coordinates
560            m[0] += xllcorner
561            m[1] += yllcorner
562           
563            # Register point and index
564            self.midpoint_coordinates[i,:] = m
565
566            # Register index of this boundary edge for use with evaluate
567            self.boundary_indices[(vol_id, edge_id)] = i
568
569
570        if verbose: log.critical('Initialise file_function')
571        self.F = file_function(filename, domain,
572                                   interpolation_points=self.midpoint_coordinates,
573                               time_thinning=time_thinning,
574                               use_cache=use_cache, 
575                               verbose=verbose)
576        self.domain = domain
577
578        # Test
579
580        # Here we'll flag indices outside the mesh as a warning
581        # as suggested by Joaquim Luis in sourceforge posting
582        # November 2007
583        # We won't make it an error as it is conceivable that
584        # only part of mesh boundary is actually used with a given
585        # file boundary sww file.
586        if (hasattr(self.F, 'indices_outside_mesh') and
587               len(self.F.indices_outside_mesh)) > 0:
588            msg = 'WARNING: File_boundary has points outside the mesh '
589            msg += 'given in %s. ' % filename
590            msg += 'See warning message issued by Interpolation_function '
591            msg += 'for details (should appear above somewhere if '
592            msg += 'verbose is True).\n'
593            msg += 'This is perfectly OK as long as the points that are '
594            msg += 'outside aren\'t used on the actual boundary segment.'
595            if verbose is True:           
596                log.critical(msg)
597            #raise Exception(msg)
598
599        q = self.F(0, point_id=0)
600
601        d = len(domain.conserved_quantities)
602        msg = 'Values specified in file %s must be ' % filename
603        msg += 'a list or an array of length %d' % d
604        assert len(q) == d, msg
605
606
607    def __repr__(self):
608        return 'File boundary'
609
610
611    def evaluate(self, vol_id=None, edge_id=None):
612        """Return linearly interpolated values based on domain.time
613               at midpoint of segment defined by vol_id and edge_id.
614        """
615
616        q = self.domain.get_conserved_quantities(vol_id, edge=edge_id)
617        t = self.domain.time
618
619        if vol_id is not None and edge_id is not None:
620            i = self.boundary_indices[vol_id, edge_id]
621            res = self.F(t, point_id=i)
622
623            if res == NAN:
624                x,y = self.midpoint_coordinates[i,:]
625                msg = 'NAN value found in file_boundary at '
626                msg += 'point id #%d: (%.2f, %.2f).\n' % (i, x, y)
627
628                if (hasattr(self.F, 'indices_outside_mesh') and
629                       len(self.F.indices_outside_mesh) > 0):
630                    # Check if NAN point is due it being outside
631                    # boundary defined in sww file.
632
633                    if i in self.F.indices_outside_mesh:
634                        msg += 'This point refers to one outside the '
635                        msg += 'mesh defined by the file %s.\n' % self.F.filename
636                        msg += 'Make sure that the file covers '
637                        msg += 'the boundary segment it is assigned to '
638                        msg += 'in set_boundary.'
639                    else:
640                        msg += 'This point is inside the mesh defined '
641                        msg += 'the file %s.\n' % self.F.filename
642                        msg += 'Check this file for NANs.'
643                raise Exception, msg
644           
645            q[0] = res[0] # Take stage, leave momentum alone
646            return q
647            #return res
648        else:
649            # raise 'Boundary call without point_id not implemented'
650            # FIXME: What should the semantics be?
651            return self.F(t)
652
653
Note: See TracBrowser for help on using the repository browser.